diff --git a/R-proj/NAMESPACE b/R-proj/NAMESPACE index f722d4e7f..d35a2450a 100644 --- a/R-proj/NAMESPACE +++ b/R-proj/NAMESPACE @@ -14,7 +14,7 @@ export(gen_rand_zonotope) export(gen_simplex) export(gen_skinny_cube) export(inner_ball) -export(rand_rotate) +export(rotate_polytope) export(round_polytope) export(sample_points) export(volume) diff --git a/R-proj/R/RcppExports.R b/R-proj/R/RcppExports.R index 98c0985f4..3dde5e511 100644 --- a/R-proj/R/RcppExports.R +++ b/R-proj/R/RcppExports.R @@ -115,42 +115,42 @@ inner_ball <- function(P) { #' #' @param kind_gen An integer to declare the type of the polytope. #' @param Vpoly_gen A boolean parameter to declare if the requested polytope has to be in V-representation. +#' @param Zono_gen A boolean parameter to declare if the requested polytope has to be a zonotope. #' @param dim_gen An integer to declare the dimension of the requested polytope. -#' @param m_gen An integer to declare the number of generators for the requested random zonotope. +#' @param m_gen An integer to declare the number of generators for the requested random zonotope or the number of vertices for a V-polytope. +#' @param seed This variable can be used to set a seed for the random polytope generator. #' #' @section warning: #' Do not use this function. #' #' @return A numerical matrix describing the requested polytope -poly_gen <- function(kind_gen, Vpoly_gen, Zono_gen, dim_gen, m_gen) { - .Call(`_volesti_poly_gen`, kind_gen, Vpoly_gen, Zono_gen, dim_gen, m_gen) +poly_gen <- function(kind_gen, Vpoly_gen, Zono_gen, dim_gen, m_gen, seed = NULL) { + .Call(`_volesti_poly_gen`, kind_gen, Vpoly_gen, Zono_gen, dim_gen, m_gen, seed) } #' An internal Rccp function for the random rotation of a convex polytope #' #' @param P A convex polytope (H-, V-polytope or a zonotope). +#' @param T Optional. A rotation matrix. #' #' @section warning: #' Do not use this function. #' #' @return A matrix that describes the rotated polytope -rotating <- function(P) { - .Call(`_volesti_rotating`, P) +rotating <- function(P, T = NULL) { + .Call(`_volesti_rotating`, P, T) } #' Internal rcpp function for the rounding of a convex polytope #' #' @param P A convex polytope (H- or V-representation or zonotope). -#' @param random_walk Optional. A string that declares the random walk. -#' @param walk_length Optional. The number of the steps for the random walk. -#' @param parameters Optional. A list for the parameters of the methods: #' #' @section warning: #' Do not use this function. #' #' @return A numerical matrix that describes the rounded polytope and contains the round value. -rounding <- function(P, random_walk = NULL, walk_length = NULL, parameters = NULL) { - .Call(`_volesti_rounding`, P, random_walk, walk_length, parameters) +rounding <- function(P) { + .Call(`_volesti_rounding`, P) } #' Sample points from a convex Polytope (H-polytope, V-polytope or a zonotope) or use direct methods for uniform sampling from the unit or the canonical or an arbitrary \eqn{d}-dimensional simplex and the boundary or the interior of a \eqn{d}-dimensional hypersphere @@ -158,22 +158,27 @@ rounding <- function(P, random_walk = NULL, walk_length = NULL, parameters = NUL #' Sample n points with uniform or multidimensional spherical gaussian -centered in an internal point- target distribution. #' The \eqn{d}-dimensional unit simplex is the set of points \eqn{\vec{x}\in \R^d}, s.t.: \eqn{\sum_i x_i\leq 1}, \eqn{x_i\geq 0}. The \eqn{d}-dimensional canonical simplex is the set of points \eqn{\vec{x}\in \R^d}, s.t.: \eqn{\sum_i x_i = 1}, \eqn{x_i\geq 0}. #' -#' @param P A convex polytope. It is an object from class (a) Hpolytope or (b) Vpolytope or (c) Zonotope. -#' @param n The number of points that the function is going to sample from the convex polytope. The default value is \eqn{100}. -#' @param distribution Optional. A string that declares the target distribution: a) \code{'uniform'} for the uniform distribution or b) \code{'gaussian'} for the multidimensional spherical distribution. The default target distribution is uniform. -#' @param random_walk Optional. 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 d) \code{'BiW'} for Billiard walk. The default walk is \code{'BiW'} for the uniform distribution or \code{'CDHR'} for the Normal distribution. -#' @param walk_length Optional. The number of the steps for the random walk. The default value is \eqn{1} for \code{'BiW'} and \eqn{\lfloor 10 + d/10\rfloor} otherwise, where \eqn{d} is the dimension that the polytope lies. -#' @param exact A boolean parameter. It should be used for the uniform sampling from the boundary or the interior of a hypersphere centered at the origin or from the unit or the canonical or an arbitrary simplex. The arbitrary simplex has to be given as a V-polytope. For the rest well known convex bodies the dimension has to be declared and the type of body as well as the radius of the hypersphere. -#' @param body A string that declares the type of the body for the exact sampling: a) \code{'unit simplex'} for the unit simplex, b) \code{'canonical simplex'} for the canonical simplex, c) \code{'hypersphere'} for the boundary of a hypersphere centered at the origin, d) \code{'ball'} for the interior of a hypersphere centered at the origin. -#' @param parameters Optional. A list for the parameters of the methods: +#' @param P A convex polytope. It is an object from class (a) Hpolytope or (b) Vpolytope or (c) Zonotope or (d) an intersection of two V-polytopes. +#' @param n The number of points that the function is going to sample from the convex polytope. +#' @param random_walk Optional. A list that declares the random walk and some related parameters as follows: +#' \itemize{ +#' \item{\code{walk} }{ A string to declare the random walk: i) \code{'CDHR'} for Coordinate Directions Hit-and-Run, ii) \code{'RDHR'} for Random Directions Hit-and-Run, iii) \code{'BaW'} for Ball Walk, iv) \code{'BiW'} for Billiard walk, v) \code{'BCDHR'} boundary sampling by keeping the extreme points of CDHR or vi) \code{'BRDHR'} boundary sampling by keeping the extreme points of RDHR. The default walk is \code{'BiW'} for the uniform distribution or \code{'CDHR'} for the Gaussian distribution.} +#' \item{\code{walk_length} }{ The number of the steps for the random walk. The default value is \eqn{5} for \code{'BiW'} and \eqn{\lfloor 10 + d/10\rfloor} otherwise.} +#' \item{\code{BaW_rad} }{ The radius for the ball walk.} +#' \item{\code{L} }{The maximum length of the billiard trajectory.} +#' } +#' @param distribution Optional. A list that declares the target density and some related parameters as follows: #' \itemize{ +#' \item{\code{density}}{A string: (a) \code{'uniform'} for the uniform distribution or b) \code{'gaussian'} for the multidimensional spherical distribution. The default target distribution is uniform.} #' \item{\code{variance} }{ The variance of the multidimensional spherical gaussian. The default value is 1.} +#' \item{\code{InnerPoint} }{ A \eqn{d}-dimensional numerical vector that defines a starting point in the interior of the polytope for the random walk and the mode of the Gaussian distribution. The default choice is the center of the Chebychev ball.} +#' } +#' @param known_body A list to request exact uniform sampling from special well known convex bodies through the following input parameters: +#' \itemize{ +#' \item{\code{body} }{ A string that declares the type of the body for the exact sampling: a) \code{'unit simplex'} for the unit simplex, b) \code{'canonical simplex'} for the canonical simplex, c) \code{'hypersphere'} for the boundary of a hypersphere centered at the origin, d) \code{'ball'} for the interior of a hypersphere centered at the origin.} #' \item{\code{dimension} }{ An integer that declares the dimension when exact sampling is enabled for a simplex or a hypersphere.} #' \item{\code{radius} }{ The radius of the \eqn{d}-dimensional hypersphere. The default value is \eqn{1}.} -#' \item{\code{BW_rad} }{ The radius for the ball walk.} -#' \item{\code{L} }{The maximum length of the billiard trajectory.} #' } -#' @param InnerPoint A \eqn{d}-dimensional numerical vector that defines a point in the interior of polytope P. #' #' @references \cite{R.Y. Rubinstein and B. Melamed, #' \dQuote{Modern simulation and modeling} \emph{ Wiley Series in Probability and Statistics,} 1998.} @@ -186,24 +191,23 @@ rounding <- function(P, random_walk = NULL, walk_length = NULL, parameters = NUL #' @examples #' # uniform distribution from the 3d unit cube in V-representation using ball walk #' P = gen_cube(3, 'V') -#' points = sample_points(P, random_walk = "BaW", walk_length = 5) +#' points = sample_points(P, n = 100, random_walk = list("walk" = "BaW", "walk_length" = 5)) #' #' # gaussian distribution from the 2d unit simplex in H-representation with variance = 2 #' A = matrix(c(-1,0,0,-1,1,1), ncol=2, nrow=3, byrow=TRUE) #' b = c(0,0,1) #' P = Hpolytope$new(A,b) -#' points = sample_points(P, distribution = "gaussian", parameters = list("variance" = 2)) +#' points = sample_points(P, n = 100, distribution = list("density" = "gaussian", "variance" = 2)) #' -#' # uniform points from the boundary of a 10-dimensional hypersphere -#' points = sample_points(exact = TRUE, body = "hypersphere", parameters = list("dimension" = 10)) +#' # uniform points from the boundary of a 2-dimensional random H-polytope +#' P = gen_rand_hpoly(2,20) +#' points = sample_points(P, n = 5000, random_walk = list("walk" = "BRDHR")) #' -#' # 10000 uniform points from a 2-d arbitrary simplex -#' V = matrix(c(2,3,-1,7,0,0),ncol = 2, nrow = 3, byrow = TRUE) -#' P = Vpolytope$new(V) -#' points = sample_points(P, n = 10000, exact = TRUE) +#' # 100 uniform points from the 2-d unit ball +#' points = sample_points(n = 100, known_body = list("body" = "ball", "dimension" = 2)) #' @export -sample_points <- function(P = NULL, n = NULL, distribution = NULL, random_walk = NULL, walk_length = NULL, exact = NULL, body = NULL, parameters = NULL, InnerPoint = NULL) { - .Call(`_volesti_sample_points`, P, n, distribution, random_walk, walk_length, exact, body, parameters, InnerPoint) +sample_points <- function(P = NULL, n = NULL, random_walk = NULL, distribution = NULL, known_body = NULL) { + .Call(`_volesti_sample_points`, P, n, random_walk, distribution, known_body) } #' The main function for volume approximation of a convex Polytope (H-polytope, V-polytope or a zonotope) @@ -211,30 +215,30 @@ sample_points <- function(P = NULL, n = NULL, distribution = NULL, random_walk = #' For the volume approximation can be used two algorithms. Either SequenceOfBalls or CoolingGaussian. A H-polytope with \eqn{m} facets is described by a \eqn{m\times d} matrix \eqn{A} and a \eqn{m}-dimensional vector \eqn{b}, s.t.: \eqn{Ax\leq b}. A V-polytope is defined as the convex hull of \eqn{m} \eqn{d}-dimensional points which correspond to the vertices of P. A zonotope is desrcibed by the Minkowski sum of \eqn{m} \eqn{d}-dimensional segments. #' #' @param P A convex polytope. It is an object from class (a) Hpolytope or (b) Vpolytope or (c) Zonotope. -#' @param walk_length Optional. The number of the steps for the random walk. The default value is \eqn{\lfloor 10 + d/10\rfloor} for SequenceOfBalls and \eqn{1} otherwise. -#' @param error Optional. Declare the upper bound for the approximation error. The default value is \eqn{1} for SequenceOfBalls and \eqn{0.1} otherwise. -#' @param inner_ball Optional. A \eqn{d+1} 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. -#' @param algo Optional. A string that declares which algorithm to use: a) \code{'SoB'} for SequenceOfBalls or b) \code{'CG'} for CoolingGaussian or c) \code{'CB'} for cooling bodies. -#' @param random_walk Optional. 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. -#' @param rounding Optional. A boolean parameter for rounding. The default value is \code{TRUE} for V-polytopes and \code{FALSE} otherwise. -#' @param parameters Optional. A list for the parameters of the algorithms: +#' @param algo Optional. A list that declares which algorithm, random walk and values of parameters to use, as follows: #' \itemize{ -#' \item{\code{Window} }{ 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{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 stopping criterion.} -#' \item{\code{L} }{The maximum length of the billiard trajectory.} +#' \item{\code{algorithm} }{ A string to set the algorithm to use: a) \code{'SoB'} for SequenceOfBalls or b) \code{'CG'} for CoolingGaussian or c) \code{'CB'} for cooling bodies. The defalut algorithm for H-polytopes is \code{'CB'} when \eqn{d\leq 200} and \code{'CG'} when \eqn{d>200}. For the other representations the default algorithm is \code{'CB'}.} +#' \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{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. #' #' @references \cite{I.Z.Emiris and V. Fisikopoulos, #' \dQuote{Practical polytope volume approximation,} \emph{ACM Trans. Math. Soft.,} 2014.}, @@ -251,32 +255,27 @@ sample_points <- function(P = NULL, n = NULL, distribution = NULL, random_walk = #' #' # calling CG algorithm for a V-polytope (3d simplex) #' P = gen_simplex(2,'V') -#' vol = volume(P, algo = "CG") +#' vol = volume(P, algo = list("algorithm" = "CG")) #' #' # calling CG algorithm for a 2-dimensional zonotope defined as the Minkowski sum of 4 segments #' Z = gen_rand_zonotope(2, 4) -#' vol = volume(Z, random_walk = "RDHR", walk_length = 5) +#' vol = volume(Z, algo = list("random_walk" = "RDHR", "walk_length" = 5)) #' @export -volume <- function(P, walk_length = NULL, error = NULL, inner_ball = NULL, algo = NULL, random_walk = NULL, rounding = NULL, parameters = NULL) { - .Call(`_volesti_volume`, P, walk_length, error, inner_ball, algo, random_walk, rounding, parameters) +volume <- function(P, algo = NULL, rounding = NULL) { + .Call(`_volesti_volume`, P, algo, rounding) } #' 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 walk_length Optional. The number of the steps for the random walk. The default value is \eqn{\lfloor 10 + d/10\rfloor} for SequenceOfBalls and \eqn{1} for CoolingGaussian. -#' @param error Optional. Declare the upper bound for the approximation error. The default value is \eqn{1} for SequenceOfBalls and \eqn{0.1} for CoolingGaussian. -#' @param inner_ball Optional. A \eqn{d+1} 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. -#' @param random_walk Optional. 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 or c) \code{'BW'} for Ball Walk. The default walk is \code{'CDHR'}. -#' @param rounding Optional. A boolean parameter for rounding. The default value is \code{FALSE}. -#' @param parameters Optional. A list for the parameters of the volume algorithm +#' @param algo_parameters Optional. A list that declares the values of the parameters of CB algorithm. #' #' @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, walk_length = NULL, error = NULL, inner_ball = NULL, random_walk = NULL, rounding = NULL, parameters = NULL) { - .Call(`_volesti_zono_approx`, Z, fit_ratio, walk_length, error, inner_ball, random_walk, rounding, parameters) +zono_approx <- function(Z, fit_ratio = NULL, algo_parameters = NULL) { + .Call(`_volesti_zono_approx`, Z, fit_ratio, algo_parameters) } diff --git a/R-proj/R/compute_indicators.R b/R-proj/R/compute_indicators.R index fd859ec18..a93199b81 100644 --- a/R-proj/R/compute_indicators.R +++ b/R-proj/R/compute_indicators.R @@ -8,6 +8,8 @@ #' @param win_length Optional. The length of the sliding window. The default value is 60. #' @param m Optional. The number of slices for the copula. The default value is 100. #' @param n Optional. The number of points to sample. The default value is \eqn{5\cdot 10^5}. +#' @param nwarning Optional. The number of consecutive indicators larger than 1 required to declare a warning period. The default value is 60. +#' @param ncrisis Optional. The number of consecutive indicators larger than 1 required to declare a crisis period. The default value is 100. #' #' @references \cite{L. Cales, A. Chalkis, I.Z. Emiris, V. Fisikopoulos, #' \dQuote{Practical volume computation of structured convex bodies, and an application to modeling portfolio dependencies and financial crises,} \emph{Proc. of Symposium on Computational Geometry, Budapest, Hungary,} 2018.} @@ -15,11 +17,13 @@ #' @return A list that contains the indicators and the corresponding vector that label each time period with respect to the market state: a) normal, b) crisis, c) warning. #' #' @export -compute_indicators <- function(returns, win_length = NULL, m = NULL, n = NULL) { +compute_indicators <- function(returns, win_length = NULL, m = NULL, n = NULL, nwarning = NULL, ncrisis = NULL) { if (is.null(win_length)) win_length = 60 if (is.null(m)) m = 100 if (is.null(n)) n = 500000 + if (is.null(nwarning)) nwarning = 60 + if (is.null(ncrisis)) ncrisis = 100 nrows = dim(returns)[1] nassets = dim(returns)[2] @@ -72,9 +76,9 @@ compute_indicators <- function(returns, win_length = NULL, m = NULL, n = NULL) { set_index = TRUE } else if (indicators[i]<1) { if(set_index){ - if(i-index+1>30 && i-index+1<60){ + if(i-index+1 > nwarning && i-index+1 < ncrisis){ col[index:i] = "warning" - } else if(i-index+1>60) { + } else if(i-index+1 > ncrisis) { col[index:i] = "crisis" } } diff --git a/R-proj/R/gen_rand_hpoly.R b/R-proj/R/gen_rand_hpoly.R index 2644529c5..354e9d268 100644 --- a/R-proj/R/gen_rand_hpoly.R +++ b/R-proj/R/gen_rand_hpoly.R @@ -4,18 +4,19 @@ #' #' @param dimension The dimension of the convex polytope. #' @param nfacets The number of the facets. +#' @param seed Optional. A fixed seed for the generator. #' #' @return A polytope class representing a H-polytope. #' @examples #' # generate a 10-dimensional polytope with 50 facets #' P = gen_rand_hpoly(10, 50) #' @export -gen_rand_hpoly <- function(dimension, nfacets) { +gen_rand_hpoly <- function(dimension, nfacets, seed = NULL) { kind_gen = 6 Vpoly_gen = FALSE - Mat = poly_gen(kind_gen, Vpoly_gen, FALSE, dimension, nfacets) + Mat = poly_gen(kind_gen, Vpoly_gen, FALSE, dimension, nfacets, seed) # first column is the vector b b = Mat[,1] diff --git a/R-proj/R/gen_rand_vpoly.R b/R-proj/R/gen_rand_vpoly.R index 9221a72cf..3711938be 100644 --- a/R-proj/R/gen_rand_vpoly.R +++ b/R-proj/R/gen_rand_vpoly.R @@ -5,13 +5,14 @@ #' @param dimension The dimension of the convex polytope. #' @param nvertices The number of the vertices. #' @param generator The body that the generator samples uniformly the vertices from: (a) 'cube' or (b) 'sphere'. +#' @param seed Optional. A fixed seed for the generator. #' #' @return A polytope class representing a V-polytope. #' @examples #' # generate a 10-dimensional polytope defined as the convex hull of 25 random vertices #' P = gen_rand_vpoly(10, 25) #' @export -gen_rand_vpoly <- function(dimension, nvertices, generator = NULL) { +gen_rand_vpoly <- function(dimension, nvertices, generator = NULL, seed = NULL) { kind_gen = 4 @@ -23,7 +24,7 @@ gen_rand_vpoly <- function(dimension, nvertices, generator = NULL) { } } - Mat = poly_gen(kind_gen, TRUE, FALSE, dimension, nvertices) + Mat = poly_gen(kind_gen, TRUE, FALSE, dimension, nvertices, seed) # first column is the vector b b = Mat[,1] diff --git a/R-proj/R/gen_rand_zonotope.R b/R-proj/R/gen_rand_zonotope.R index 399ccdccf..d4ea92eb2 100644 --- a/R-proj/R/gen_rand_zonotope.R +++ b/R-proj/R/gen_rand_zonotope.R @@ -5,14 +5,15 @@ #' @param dimension The dimension of the zonotope. #' @param nsegments The number of segments that generate the zonotope. #' @param generator The distribution to pick the length of each segment from \eqn{[0,100]}: (a) 'uniform', (b) 'gaussian' or (c) 'exponential'. -#' +#' @param seed Optional. A fixed seed for the generator. +#' #' @return A polytope class representing a zonotope. #' #' @examples #' # generate a 10-dimensional zonotope defined by the Minkowski sum of 20 segments #' P = gen_rand_zonotope(10, 20) #' @export -gen_rand_zonotope <- function(dimension, nsegments, generator = NULL) { +gen_rand_zonotope <- function(dimension, nsegments, generator = NULL, seed = NULL) { kind_gen = 1 @@ -26,7 +27,7 @@ gen_rand_zonotope <- function(dimension, nsegments, generator = NULL) { } } - Mat = poly_gen(kind_gen, FALSE, TRUE, dimension, nsegments) + Mat = poly_gen(kind_gen, FALSE, TRUE, dimension, nsegments, seed) # first column is the vector b b = Mat[,1] diff --git a/R-proj/R/rand_rotate.R b/R-proj/R/rotate_polytope.R similarity index 73% rename from R-proj/R/rand_rotate.R rename to R-proj/R/rotate_polytope.R index fae4d8841..3a20276c0 100644 --- a/R-proj/R/rand_rotate.R +++ b/R-proj/R/rotate_polytope.R @@ -1,10 +1,11 @@ -#' Apply a random rotation to a convex polytope (H-polytope, V-polytope or a zonotope) +#' Apply a random rotation to a convex polytope (H-polytope, V-polytope, zonotope or intersection of two V-polytopes) #' -#' Given a convex H or V polytope or a zonotope as input this function applies a random rotation. +#' Given a convex H- or V- polytope or a zonotope or an intersection of two V-polytopes as input, this function applies (a) a random rotation or (b) a given rotation by an input matrix \eqn{T}. #' -#' @param P A convex polytope. It is an object from class (a) Hpolytope or (b) Vpolytope or (c) Zonotope. +#' @param P A convex polytope. It is an object from class (a) Hpolytope, (b) Vpolytope, (c) Zonotope, (d) intersection of two V-polytopes. +#' @param T Optional. A \eqn{d\times d} rotation matrix. #' -#' @return A list that contains the rotated polytope and the matrix of the linear transformation. +#' @return A list that contains the rotated polytope and the matrix \eqn{T} of the linear transformation. #' #' @details Let \eqn{P} be the given polytope and \eqn{Q} the rotated one and \eqn{T} be the matrix of the linear transformation. #' \itemize{ @@ -16,20 +17,20 @@ #' @examples #' # rotate a H-polytope (2d unit simplex) #' P = gen_simplex(2,'H') -#' poly_matrix_list = rand_rotate(P) +#' poly_matrix_list = rotate_polytope(P) #' #' # rotate a V-polytope (3d cube) #' P = gen_cube(3, 'V') -#' poly_matrix_list = rand_rotate(P) +#' poly_matrix_list = rotate_polytope(P) #' #' # rotate a 5-dimensional zonotope defined by the Minkowski sum of 15 segments #' Z = gen_rand_zonotope(3,6) -#' poly_matrix_list = rand_rotate(Z) +#' poly_matrix_list = rotate_polytope(Z) #' @export -rand_rotate <- function(P){ +rotate_polytope <- function(P, T = NULL){ #call rcpp rotating function - Mat = rotating(P) + Mat = rotating(P, T) n = P$dimension m=dim(Mat)[2]-n diff --git a/R-proj/R/round_polytope.R b/R-proj/R/round_polytope.R index 16b8367ae..524aab22b 100644 --- a/R-proj/R/round_polytope.R +++ b/R-proj/R/round_polytope.R @@ -3,16 +3,12 @@ #' 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 random_walk Optional. 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 d) \code{'BiW'} for Billiard walk. The default walk is \code{'CDHR'} for H-polytope and \code{'BiW'} for the other representations. -#' @param walk_length Optional. The number of the steps for the random walk. The default value is \eqn{1} for \code{'BiW'} and \eqn{\lfloor 10 + d/10\rfloor} otherwise. -#' @param parameters Optional. A list for the parameters of the methods: -#' \itemize{ -#' \item{\code{BW_rad} }{ The radius for the ball walk.} -#' \item{\code{diameter} }{ The diameter of the polytope. It is used to set the maximum length of the trajectory in billiard walk.} -#' } #' #' @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. #' +#' @references \cite{Michael J. Todd and E. Alper Yildirim, +#' \dQuote{On Khachiyan’s Algorithm for the Computation of Minimum Volume Enclosing Ellipsoids,} \emph{Discrete Applied Mathematics,} 2007.} +#' #' @examples #' # rotate a H-polytope (2d unit simplex) #' A = matrix(c(-1,0,0,-1,1,1), ncol=2, nrow=3, byrow=TRUE) @@ -22,15 +18,15 @@ #' #' # rotate a V-polytope (3d unit cube) using Random Directions HnR with step equal to 50 #' P = gen_cube(3, 'V') -#' ListVpoly = round_polytope(P, random_walk = 'RDHR', walk_length = 50) +#' ListVpoly = round_polytope(P) #' #' # round a 2-dimensional zonotope defined by 6 generators using ball walk #' Z = gen_rand_zonotope(2,6) -#' ListZono = round_polytope(Z, random_walk = 'BW') +#' ListZono = round_polytope(Z) #' @export -round_polytope <- function(P, random_walk = NULL, walk_length = NULL, parameters = NULL){ +round_polytope <- function(P){ - ret_list = rounding(P, random_walk, walk_length, parameters) + ret_list = rounding(P) #get the matrix that describes the polytope Mat = ret_list$Mat diff --git a/R-proj/R/zonotope_approximation.R b/R-proj/R/zonotope_approximation.R index d8b38a1e8..e17ca94e2 100644 --- a/R-proj/R/zonotope_approximation.R +++ b/R-proj/R/zonotope_approximation.R @@ -4,39 +4,37 @@ #' #' @param Z A zonotope. #' @param fit_ratio Optional. A boolean parameter to request the computation of the ratio of fitness. -#' @param walk_length Optional. The number of the steps for the random walk. The default value is \eqn{\lfloor 10 + d/10\rfloor} for SequenceOfBalls and \eqn{1} for CoolingGaussian. -#' @param error Optional. Declare the upper bound for the approximation error. The default value is \eqn{1} for SequenceOfBalls and \eqn{0.1} for CoolingGaussian. -#' @param inner_ball Optional. A \eqn{d+1} 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. -#' @param random_walk Optional. 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 or c) \code{'BW'} for Ball Walk. The default walk is \code{'CDHR'}. -#' @param rounding Optional. A boolean parameter for rounding. The default value is \code{FALSE}. -#' @param parameters Optional. A list for the parameters of the algorithms: +#' @param algo_parameters Optional. A list that declares the values of the parameters of CB algorithm as follows: #' \itemize{ -#' \item{\code{Window} }{ 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 BAN algorithm. The default value is \eqn{0.1}.} -#' \item{\code{lb} }{ The upper bound for the ratios in MMC in BAN 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.} -#' \item{\code{nu} }{ The degrees of freedom for the t-student distribution in t-tests in BAN algorithm. The default value is \eqn{10}.} -#' \item{\code{alpha} }{ The significance level for the t-tests in BAN algorithm. The default values is 0.2.} -#' \item{\code{prob} }{ The probability is used for the empirical confidence interval in ratio estimation of BAN algorithm. The default value is \eqn{0.75}.} -#' \item{\code{hpoly} }{ A boolean parameter to use H-polytopes in MMC of BAN algorithm. The default value is \code{FALSE}.} -#' \item{\code{minmaxW} }{ A boolean parameter to use the sliding window with a minmax values stopping criterion.} -#' \item{\code{L} }{The maximum length of the billiard trajectory.} +#' \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{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.} #' } #' #' @return A list that contains the approximation body in H-representation and the ratio of fitness #' #' @examples #' # over-approximate a 2-dimensional zonotope with 10 generators and compute the ratio of fitness -#' Z = GenZonotope(2,10) +#' Z = gen_rand_zonotope(2,10) #' retList = zonotope_approximation(Z = Z, fit_ratio = TRUE) #' #' @export -zonotope_approximation <- function(Z, fit_ratio = NULL, walk_length = NULL, error = NULL, - inner_ball = NULL, random_walk = NULL, rounding = NULL, - parameters = NULL){ +zonotope_approximation <- function(Z, fit_ratio = NULL, algo_parameters = NULL){ - ret_list = zono_approx(Z, fit_ratio, walk_length, error, inner_ball, random_walk, rounding, parameters) + ret_list = zono_approx(Z, fit_ratio, algo_parameters) Mat = ret_list$Mat diff --git a/R-proj/man/compute_indicators.Rd b/R-proj/man/compute_indicators.Rd index bf9d720b4..df93bc806 100644 --- a/R-proj/man/compute_indicators.Rd +++ b/R-proj/man/compute_indicators.Rd @@ -4,7 +4,8 @@ \alias{compute_indicators} \title{Compute an indicator for each time period that describes the state of a market.} \usage{ -compute_indicators(returns, win_length = NULL, m = NULL, n = NULL) +compute_indicators(returns, win_length = NULL, m = NULL, n = NULL, + nwarning = NULL, ncrisis = NULL) } \arguments{ \item{returns}{A \eqn{d}-dimensional vector that describes the direction of the first family of parallel hyperplanes.} @@ -14,6 +15,10 @@ compute_indicators(returns, win_length = NULL, m = NULL, n = NULL) \item{m}{Optional. The number of slices for the copula. The default value is 100.} \item{n}{Optional. The number of points to sample. The default value is \eqn{5\cdot 10^5}.} + +\item{nwarning}{Optional. The number of consecutive indicators larger than 1 required to declare a warning period. The default value is 60.} + +\item{ncrisis}{Optional. The number of consecutive indicators larger than 1 required to declare a crisis period. The default value is 100.} } \value{ A list that contains the indicators and the corresponding vector that label each time period with respect to the market state: a) normal, b) crisis, c) warning. diff --git a/R-proj/man/gen_rand_hpoly.Rd b/R-proj/man/gen_rand_hpoly.Rd index 0f72d562a..b5c9adf39 100644 --- a/R-proj/man/gen_rand_hpoly.Rd +++ b/R-proj/man/gen_rand_hpoly.Rd @@ -4,12 +4,14 @@ \alias{gen_rand_hpoly} \title{Generator function for random H-polytopes} \usage{ -gen_rand_hpoly(dimension, nfacets) +gen_rand_hpoly(dimension, nfacets, seed = NULL) } \arguments{ \item{dimension}{The dimension of the convex polytope.} \item{nfacets}{The number of the facets.} + +\item{seed}{Optional. A fixed seed for the generator.} } \value{ A polytope class representing a H-polytope. diff --git a/R-proj/man/gen_rand_vpoly.Rd b/R-proj/man/gen_rand_vpoly.Rd index d57b56c33..95dffb981 100644 --- a/R-proj/man/gen_rand_vpoly.Rd +++ b/R-proj/man/gen_rand_vpoly.Rd @@ -4,7 +4,7 @@ \alias{gen_rand_vpoly} \title{Generator function for random V-polytopes} \usage{ -gen_rand_vpoly(dimension, nvertices, generator = NULL) +gen_rand_vpoly(dimension, nvertices, generator = NULL, seed = NULL) } \arguments{ \item{dimension}{The dimension of the convex polytope.} @@ -12,6 +12,8 @@ gen_rand_vpoly(dimension, nvertices, generator = NULL) \item{nvertices}{The number of the vertices.} \item{generator}{The body that the generator samples uniformly the vertices from: (a) 'cube' or (b) 'sphere'.} + +\item{seed}{Optional. A fixed seed for the generator.} } \value{ A polytope class representing a V-polytope. diff --git a/R-proj/man/gen_rand_zonotope.Rd b/R-proj/man/gen_rand_zonotope.Rd index dca32115b..bbe3e19cc 100644 --- a/R-proj/man/gen_rand_zonotope.Rd +++ b/R-proj/man/gen_rand_zonotope.Rd @@ -4,7 +4,7 @@ \alias{gen_rand_zonotope} \title{Generator function for zonotopes} \usage{ -gen_rand_zonotope(dimension, nsegments, generator = NULL) +gen_rand_zonotope(dimension, nsegments, generator = NULL, seed = NULL) } \arguments{ \item{dimension}{The dimension of the zonotope.} @@ -12,6 +12,8 @@ gen_rand_zonotope(dimension, nsegments, generator = NULL) \item{nsegments}{The number of segments that generate the zonotope.} \item{generator}{The distribution to pick the length of each segment from \eqn{[0,100]}: (a) 'uniform', (b) 'gaussian' or (c) 'exponential'.} + +\item{seed}{Optional. A fixed seed for the generator.} } \value{ A polytope class representing a zonotope. diff --git a/R-proj/man/poly_gen.Rd b/R-proj/man/poly_gen.Rd index 2e7ce5e02..e16aae1ee 100644 --- a/R-proj/man/poly_gen.Rd +++ b/R-proj/man/poly_gen.Rd @@ -4,16 +4,20 @@ \alias{poly_gen} \title{An internal Rccp function as a polytope generator} \usage{ -poly_gen(kind_gen, Vpoly_gen, Zono_gen, dim_gen, m_gen) +poly_gen(kind_gen, Vpoly_gen, Zono_gen, dim_gen, m_gen, seed = NULL) } \arguments{ \item{kind_gen}{An integer to declare the type of the polytope.} \item{Vpoly_gen}{A boolean parameter to declare if the requested polytope has to be in V-representation.} +\item{Zono_gen}{A boolean parameter to declare if the requested polytope has to be a zonotope.} + \item{dim_gen}{An integer to declare the dimension of the requested polytope.} -\item{m_gen}{An integer to declare the number of generators for the requested random zonotope.} +\item{m_gen}{An integer to declare the number of generators for the requested random zonotope or the number of vertices for a V-polytope.} + +\item{seed}{This variable can be used to set a seed for the random polytope generator.} } \value{ A numerical matrix describing the requested polytope diff --git a/R-proj/man/rand_rotate.Rd b/R-proj/man/rotate_polytope.Rd similarity index 63% rename from R-proj/man/rand_rotate.Rd rename to R-proj/man/rotate_polytope.Rd index 574953868..c53583536 100644 --- a/R-proj/man/rand_rotate.Rd +++ b/R-proj/man/rotate_polytope.Rd @@ -1,19 +1,21 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/rand_rotate.R -\name{rand_rotate} -\alias{rand_rotate} -\title{Apply a random rotation to a convex polytope (H-polytope, V-polytope or a zonotope)} +% Please edit documentation in R/rotate_polytope.R +\name{rotate_polytope} +\alias{rotate_polytope} +\title{Apply a random rotation to a convex polytope (H-polytope, V-polytope, zonotope or intersection of two V-polytopes)} \usage{ -rand_rotate(P) +rotate_polytope(P, T = NULL) } \arguments{ -\item{P}{A convex polytope. It is an object from class (a) Hpolytope or (b) Vpolytope or (c) Zonotope.} +\item{P}{A convex polytope. It is an object from class (a) Hpolytope, (b) Vpolytope, (c) Zonotope, (d) intersection of two V-polytopes.} + +\item{T}{Optional. A \eqn{d\times d} rotation matrix.} } \value{ -A list that contains the rotated polytope and the matrix of the linear transformation. +A list that contains the rotated polytope and the matrix \eqn{T} of the linear transformation. } \description{ -Given a convex H or V polytope or a zonotope as input this function applies a random rotation. +Given a convex H- or V- polytope or a zonotope or an intersection of two V-polytopes as input, this function applies (a) a random rotation or (b) a given rotation by an input matrix \eqn{T}. } \details{ Let \eqn{P} be the given polytope and \eqn{Q} the rotated one and \eqn{T} be the matrix of the linear transformation. @@ -27,13 +29,13 @@ Let \eqn{P} be the given polytope and \eqn{Q} the rotated one and \eqn{T} be the \examples{ # rotate a H-polytope (2d unit simplex) P = gen_simplex(2,'H') -poly_matrix_list = rand_rotate(P) +poly_matrix_list = rotate_polytope(P) # rotate a V-polytope (3d cube) P = gen_cube(3, 'V') -poly_matrix_list = rand_rotate(P) +poly_matrix_list = rotate_polytope(P) # rotate a 5-dimensional zonotope defined by the Minkowski sum of 15 segments Z = gen_rand_zonotope(3,6) -poly_matrix_list = rand_rotate(Z) +poly_matrix_list = rotate_polytope(Z) } diff --git a/R-proj/man/rotating.Rd b/R-proj/man/rotating.Rd index dc99d9f2d..f3514fb09 100644 --- a/R-proj/man/rotating.Rd +++ b/R-proj/man/rotating.Rd @@ -4,10 +4,12 @@ \alias{rotating} \title{An internal Rccp function for the random rotation of a convex polytope} \usage{ -rotating(P) +rotating(P, T = NULL) } \arguments{ \item{P}{A convex polytope (H-, V-polytope or a zonotope).} + +\item{T}{Optional. A rotation matrix.} } \value{ A matrix that describes the rotated polytope diff --git a/R-proj/man/round_polytope.Rd b/R-proj/man/round_polytope.Rd index 627db872a..6874f8eff 100644 --- a/R-proj/man/round_polytope.Rd +++ b/R-proj/man/round_polytope.Rd @@ -4,21 +4,10 @@ \alias{round_polytope} \title{Apply rounding to a convex polytope (H-polytope, V-polytope or a zonotope)} \usage{ -round_polytope(P, random_walk = NULL, walk_length = NULL, - parameters = NULL) +round_polytope(P) } \arguments{ \item{P}{A convex polytope. It is an object from class (a) Hpolytope or (b) Vpolytope or (c) Zonotope.} - -\item{random_walk}{Optional. 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 d) \code{'BiW'} for Billiard walk. The default walk is \code{'CDHR'} for H-polytope and \code{'BiW'} for the other representations.} - -\item{walk_length}{Optional. The number of the steps for the random walk. The default value is \eqn{1} for \code{'BiW'} and \eqn{\lfloor 10 + d/10\rfloor} otherwise.} - -\item{parameters}{Optional. A list for the parameters of the methods: -\itemize{ - \item{\code{BW_rad} }{ The radius for the ball walk.} - \item{\code{diameter} }{ The diameter of the polytope. It is used to set the maximum length of the trajectory in billiard walk.} -}} } \value{ 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. @@ -35,9 +24,13 @@ listHpoly = round_polytope(P) # rotate a V-polytope (3d unit cube) using Random Directions HnR with step equal to 50 P = gen_cube(3, 'V') -ListVpoly = round_polytope(P, random_walk = 'RDHR', walk_length = 50) +ListVpoly = round_polytope(P) # round a 2-dimensional zonotope defined by 6 generators using ball walk Z = gen_rand_zonotope(2,6) -ListZono = round_polytope(Z, random_walk = 'BW') +ListZono = round_polytope(Z) +} +\references{ +\cite{Michael J. Todd and E. Alper Yildirim, +\dQuote{On Khachiyan’s Algorithm for the Computation of Minimum Volume Enclosing Ellipsoids,} \emph{Discrete Applied Mathematics,} 2007.} } diff --git a/R-proj/man/rounding.Rd b/R-proj/man/rounding.Rd index 10e5a6564..e28fdc0e9 100644 --- a/R-proj/man/rounding.Rd +++ b/R-proj/man/rounding.Rd @@ -4,17 +4,10 @@ \alias{rounding} \title{Internal rcpp function for the rounding of a convex polytope} \usage{ -rounding(P, random_walk = NULL, walk_length = NULL, - parameters = NULL) +rounding(P) } \arguments{ \item{P}{A convex polytope (H- or V-representation or zonotope).} - -\item{random_walk}{Optional. A string that declares the random walk.} - -\item{walk_length}{Optional. The number of the steps for the random walk.} - -\item{parameters}{Optional. A list for the parameters of the methods:} } \value{ A numerical matrix that describes the rounded polytope and contains the round value. diff --git a/R-proj/man/sample_points.Rd b/R-proj/man/sample_points.Rd index 3d18a43a9..df0f866f8 100644 --- a/R-proj/man/sample_points.Rd +++ b/R-proj/man/sample_points.Rd @@ -4,35 +4,35 @@ \alias{sample_points} \title{Sample points from a convex Polytope (H-polytope, V-polytope or a zonotope) or use direct methods for uniform sampling from the unit or the canonical or an arbitrary \eqn{d}-dimensional simplex and the boundary or the interior of a \eqn{d}-dimensional hypersphere} \usage{ -sample_points(P = NULL, n = NULL, distribution = NULL, - random_walk = NULL, walk_length = NULL, exact = NULL, - body = NULL, parameters = NULL, InnerPoint = NULL) +sample_points(P = NULL, n = NULL, random_walk = NULL, + distribution = NULL, known_body = NULL) } \arguments{ -\item{P}{A convex polytope. It is an object from class (a) Hpolytope or (b) Vpolytope or (c) Zonotope.} +\item{P}{A convex polytope. It is an object from class (a) Hpolytope or (b) Vpolytope or (c) Zonotope or (d) an intersection of two V-polytopes.} -\item{n}{The number of points that the function is going to sample from the convex polytope. The default value is \eqn{100}.} +\item{n}{The number of points that the function is going to sample from the convex polytope.} -\item{distribution}{Optional. A string that declares the target distribution: a) \code{'uniform'} for the uniform distribution or b) \code{'gaussian'} for the multidimensional spherical distribution. The default target distribution is uniform.} - -\item{random_walk}{Optional. 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 d) \code{'BiW'} for Billiard walk. The default walk is \code{'BiW'} for the uniform distribution or \code{'CDHR'} for the Normal distribution.} - -\item{walk_length}{Optional. The number of the steps for the random walk. The default value is \eqn{1} for \code{'BiW'} and \eqn{\lfloor 10 + d/10\rfloor} otherwise, where \eqn{d} is the dimension that the polytope lies.} - -\item{exact}{A boolean parameter. It should be used for the uniform sampling from the boundary or the interior of a hypersphere centered at the origin or from the unit or the canonical or an arbitrary simplex. The arbitrary simplex has to be given as a V-polytope. For the rest well known convex bodies the dimension has to be declared and the type of body as well as the radius of the hypersphere.} - -\item{body}{A string that declares the type of the body for the exact sampling: a) \code{'unit simplex'} for the unit simplex, b) \code{'canonical simplex'} for the canonical simplex, c) \code{'hypersphere'} for the boundary of a hypersphere centered at the origin, d) \code{'ball'} for the interior of a hypersphere centered at the origin.} +\item{random_walk}{Optional. A list that declares the random walk and some related parameters as follows: +\itemize{ +\item{\code{walk} }{ A string to declare the random walk: i) \code{'CDHR'} for Coordinate Directions Hit-and-Run, ii) \code{'RDHR'} for Random Directions Hit-and-Run, iii) \code{'BaW'} for Ball Walk, iv) \code{'BiW'} for Billiard walk, v) \code{'BCDHR'} boundary sampling by keeping the extreme points of CDHR or vi) \code{'BRDHR'} boundary sampling by keeping the extreme points of RDHR. The default walk is \code{'BiW'} for the uniform distribution or \code{'CDHR'} for the Gaussian distribution.} +\item{\code{walk_length} }{ The number of the steps for the random walk. The default value is \eqn{5} for \code{'BiW'} and \eqn{\lfloor 10 + d/10\rfloor} otherwise.} +\item{\code{BaW_rad} }{ The radius for the ball walk.} +\item{\code{L} }{The maximum length of the billiard trajectory.} +}} -\item{parameters}{Optional. A list for the parameters of the methods: +\item{distribution}{Optional. A list that declares the target density and some related parameters as follows: \itemize{ +\item{\code{density}}{A string: (a) \code{'uniform'} for the uniform distribution or b) \code{'gaussian'} for the multidimensional spherical distribution. The default target distribution is uniform.} \item{\code{variance} }{ The variance of the multidimensional spherical gaussian. The default value is 1.} +\item{\code{InnerPoint} }{ A \eqn{d}-dimensional numerical vector that defines a starting point in the interior of the polytope for the random walk and the mode of the Gaussian distribution. The default choice is the center of the Chebychev ball.} +}} + +\item{known_body}{A list to request exact uniform sampling from special well known convex bodies through the following input parameters: +\itemize{ +\item{\code{body} }{ A string that declares the type of the body for the exact sampling: a) \code{'unit simplex'} for the unit simplex, b) \code{'canonical simplex'} for the canonical simplex, c) \code{'hypersphere'} for the boundary of a hypersphere centered at the origin, d) \code{'ball'} for the interior of a hypersphere centered at the origin.} \item{\code{dimension} }{ An integer that declares the dimension when exact sampling is enabled for a simplex or a hypersphere.} \item{\code{radius} }{ The radius of the \eqn{d}-dimensional hypersphere. The default value is \eqn{1}.} -\item{\code{BW_rad} }{ The radius for the ball walk.} -\item{\code{L} }{The maximum length of the billiard trajectory.} }} - -\item{InnerPoint}{A \eqn{d}-dimensional numerical vector that defines a point in the interior of polytope P.} } \value{ A \eqn{d\times n} matrix that contains, column-wise, the sampled points from the convex polytope P. @@ -44,21 +44,20 @@ The \eqn{d}-dimensional unit simplex is the set of points \eqn{\vec{x}\in \R^d}, \examples{ # uniform distribution from the 3d unit cube in V-representation using ball walk P = gen_cube(3, 'V') -points = sample_points(P, random_walk = "BaW", walk_length = 5) +points = sample_points(P, n = 100, random_walk = list("walk" = "BaW", "walk_length" = 5)) # gaussian distribution from the 2d unit simplex in H-representation with variance = 2 A = matrix(c(-1,0,0,-1,1,1), ncol=2, nrow=3, byrow=TRUE) b = c(0,0,1) P = Hpolytope$new(A,b) -points = sample_points(P, distribution = "gaussian", parameters = list("variance" = 2)) +points = sample_points(P, n = 100, distribution = list("density" = "gaussian", "variance" = 2)) -# uniform points from the boundary of a 10-dimensional hypersphere -points = sample_points(exact = TRUE, body = "hypersphere", parameters = list("dimension" = 10)) +# uniform points from the boundary of a 2-dimensional random H-polytope +P = gen_rand_hpoly(2,20) +points = sample_points(P, n = 5000, random_walk = list("walk" = "BRDHR")) -# 10000 uniform points from a 2-d arbitrary simplex -V = matrix(c(2,3,-1,7,0,0),ncol = 2, nrow = 3, byrow = TRUE) -P = Vpolytope$new(V) -points = sample_points(P, n = 10000, exact = TRUE) +# 100 uniform points from the 2-d unit ball +points = sample_points(n = 100, known_body = list("body" = "ball", "dimension" = 2)) } \references{ \cite{R.Y. Rubinstein and B. Melamed, diff --git a/R-proj/man/volume.Rd b/R-proj/man/volume.Rd index f1f2b7eac..321a1ae5c 100644 --- a/R-proj/man/volume.Rd +++ b/R-proj/man/volume.Rd @@ -4,43 +4,36 @@ \alias{volume} \title{The main function for volume approximation of a convex Polytope (H-polytope, V-polytope or a zonotope)} \usage{ -volume(P, walk_length = NULL, error = NULL, inner_ball = NULL, - algo = NULL, random_walk = NULL, rounding = NULL, - parameters = NULL) +volume(P, algo = NULL, rounding = NULL) } \arguments{ \item{P}{A convex polytope. It is an object from class (a) Hpolytope or (b) Vpolytope or (c) Zonotope.} -\item{walk_length}{Optional. The number of the steps for the random walk. The default value is \eqn{\lfloor 10 + d/10\rfloor} for SequenceOfBalls and \eqn{1} otherwise.} - -\item{error}{Optional. Declare the upper bound for the approximation error. The default value is \eqn{1} for SequenceOfBalls and \eqn{0.1} otherwise.} - -\item{inner_ball}{Optional. A \eqn{d+1} 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{algo}{Optional. A string that declares which algorithm to use: a) \code{'SoB'} for SequenceOfBalls or b) \code{'CG'} for CoolingGaussian or c) \code{'CB'} for cooling bodies.} - -\item{random_walk}{Optional. 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{rounding}{Optional. A boolean parameter for rounding. The default value is \code{TRUE} for V-polytopes and \code{FALSE} otherwise.} - -\item{parameters}{Optional. A list for the parameters of the algorithms: +\item{algo}{Optional. A list that declares which algorithm, random walk and values of parameters to use, as follows: \itemize{ -\item{\code{Window} }{ 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{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 stopping criterion.} - \item{\code{L} }{The maximum length of the billiard trajectory.} +\item{\code{algorithm} }{ A string to set the algorithm to use: a) \code{'SoB'} for SequenceOfBalls or b) \code{'CG'} for CoolingGaussian or c) \code{'CB'} for cooling bodies. The defalut algorithm for H-polytopes is \code{'CB'} when \eqn{d\leq 200} and \code{'CG'} when \eqn{d>200}. For the other representations the default algorithm is \code{'CB'}.} +\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{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.} } \value{ The approximation of the volume of a convex polytope. @@ -55,11 +48,11 @@ vol = volume(P) # calling CG algorithm for a V-polytope (3d simplex) P = gen_simplex(2,'V') -vol = volume(P, algo = "CG") +vol = volume(P, algo = list("algorithm" = "CG")) # calling CG algorithm for a 2-dimensional zonotope defined as the Minkowski sum of 4 segments Z = gen_rand_zonotope(2, 4) -vol = volume(Z, random_walk = "RDHR", walk_length = 5) +vol = volume(Z, algo = list("random_walk" = "RDHR", "walk_length" = 5)) } \references{ \cite{I.Z.Emiris and V. Fisikopoulos, diff --git a/R-proj/man/zono_approx.Rd b/R-proj/man/zono_approx.Rd index 09daa5459..0123f1519 100644 --- a/R-proj/man/zono_approx.Rd +++ b/R-proj/man/zono_approx.Rd @@ -4,26 +4,14 @@ \alias{zono_approx} \title{An internal Rccp function for the over-approximation of a zonotope} \usage{ -zono_approx(Z, fit_ratio = NULL, walk_length = NULL, error = NULL, - inner_ball = NULL, random_walk = NULL, rounding = NULL, - parameters = NULL) +zono_approx(Z, fit_ratio = NULL, algo_parameters = NULL) } \arguments{ \item{Z}{A zonotope.} \item{fit_ratio}{Optional. A boolean parameter to request the computation of the ratio of fitness.} -\item{walk_length}{Optional. The number of the steps for the random walk. The default value is \eqn{\lfloor 10 + d/10\rfloor} for SequenceOfBalls and \eqn{1} for CoolingGaussian.} - -\item{error}{Optional. Declare the upper bound for the approximation error. The default value is \eqn{1} for SequenceOfBalls and \eqn{0.1} for CoolingGaussian.} - -\item{inner_ball}{Optional. A \eqn{d+1} 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{random_walk}{Optional. 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 or c) \code{'BW'} for Ball Walk. The default walk is \code{'CDHR'}.} - -\item{rounding}{Optional. A boolean parameter for rounding. The default value is \code{FALSE}.} - -\item{parameters}{Optional. A list for the parameters of the volume algorithm} +\item{algo_parameters}{Optional. A list that declares the values of the parameters of CB algorithm.} } \value{ A List that contains a numerical matrix that describes the PCA approximation as a H-polytope and the ratio of fitness. diff --git a/R-proj/man/zonotope_approximation.Rd b/R-proj/man/zonotope_approximation.Rd index 40263c74b..859155a7b 100644 --- a/R-proj/man/zonotope_approximation.Rd +++ b/R-proj/man/zonotope_approximation.Rd @@ -4,38 +4,31 @@ \alias{zonotope_approximation} \title{A function to over-approximate a zonotope with PCA method and to evaluate the approximation by computing a ratio of fitness.} \usage{ -zonotope_approximation(Z, fit_ratio = NULL, walk_length = NULL, - error = NULL, inner_ball = NULL, random_walk = NULL, - rounding = NULL, parameters = NULL) +zonotope_approximation(Z, fit_ratio = NULL, algo_parameters = NULL) } \arguments{ \item{Z}{A zonotope.} \item{fit_ratio}{Optional. A boolean parameter to request the computation of the ratio of fitness.} -\item{walk_length}{Optional. The number of the steps for the random walk. The default value is \eqn{\lfloor 10 + d/10\rfloor} for SequenceOfBalls and \eqn{1} for CoolingGaussian.} - -\item{error}{Optional. Declare the upper bound for the approximation error. The default value is \eqn{1} for SequenceOfBalls and \eqn{0.1} for CoolingGaussian.} - -\item{inner_ball}{Optional. A \eqn{d+1} 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{random_walk}{Optional. 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 or c) \code{'BW'} for Ball Walk. The default walk is \code{'CDHR'}.} - -\item{rounding}{Optional. A boolean parameter for rounding. The default value is \code{FALSE}.} - -\item{parameters}{Optional. A list for the parameters of the algorithms: +\item{algo_parameters}{Optional. A list that declares the values of the parameters of CB algorithm as follows: \itemize{ - \item{\code{Window} }{ 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 BAN algorithm. The default value is \eqn{0.1}.} - \item{\code{lb} }{ The upper bound for the ratios in MMC in BAN 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.} - \item{\code{nu} }{ The degrees of freedom for the t-student distribution in t-tests in BAN algorithm. The default value is \eqn{10}.} - \item{\code{alpha} }{ The significance level for the t-tests in BAN algorithm. The default values is 0.2.} - \item{\code{prob} }{ The probability is used for the empirical confidence interval in ratio estimation of BAN algorithm. The default value is \eqn{0.75}.} - \item{\code{hpoly} }{ A boolean parameter to use H-polytopes in MMC of BAN algorithm. The default value is \code{FALSE}.} - \item{\code{minmaxW} }{ A boolean parameter to use the sliding window with a minmax values stopping criterion.} - \item{\code{L} }{The maximum length of the billiard trajectory.} +\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{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.} }} } \value{ @@ -46,7 +39,7 @@ For the evaluation of the PCA method the exact volume of the approximation body } \examples{ # over-approximate a 2-dimensional zonotope with 10 generators and compute the ratio of fitness -Z = GenZonotope(2,10) +Z = gen_rand_zonotope(2,10) retList = zonotope_approximation(Z = Z, fit_ratio = TRUE) } diff --git a/R-proj/src/poly_gen.cpp b/R-proj/src/poly_gen.cpp index 03cea06fb..6268abe3c 100644 --- a/R-proj/src/poly_gen.cpp +++ b/R-proj/src/poly_gen.cpp @@ -28,15 +28,18 @@ //' //' @param kind_gen An integer to declare the type of the polytope. //' @param Vpoly_gen A boolean parameter to declare if the requested polytope has to be in V-representation. +//' @param Zono_gen A boolean parameter to declare if the requested polytope has to be a zonotope. //' @param dim_gen An integer to declare the dimension of the requested polytope. -//' @param m_gen An integer to declare the number of generators for the requested random zonotope. +//' @param m_gen An integer to declare the number of generators for the requested random zonotope or the number of vertices for a V-polytope. +//' @param seed This variable can be used to set a seed for the random polytope generator. //' //' @section warning: //' Do not use this function. //' //' @return A numerical matrix describing the requested polytope // [[Rcpp::export]] -Rcpp::NumericMatrix poly_gen (int kind_gen, bool Vpoly_gen, bool Zono_gen, int dim_gen, int m_gen) { +Rcpp::NumericMatrix poly_gen (int kind_gen, bool Vpoly_gen, bool Zono_gen, int dim_gen, int m_gen, + Rcpp::Nullable seed = R_NilValue) { typedef double NT; typedef Cartesian Kernel; @@ -46,16 +49,18 @@ Rcpp::NumericMatrix poly_gen (int kind_gen, bool Vpoly_gen, bool Zono_gen, int d typedef VPolytope Vpolytope; typedef Zonotope zonotope; + double seed2 = (!seed.isNotNull()) ? std::numeric_limits::signaling_NaN() : Rcpp::as(seed); + if (Zono_gen) { switch (kind_gen) { case 1: - return extractMatPoly(gen_zonotope_uniform(dim_gen, m_gen)); + return extractMatPoly(gen_zonotope_uniform(dim_gen, m_gen, seed2)); case 2: - return extractMatPoly(gen_zonotope_gaussian(dim_gen, m_gen)); + return extractMatPoly(gen_zonotope_gaussian(dim_gen, m_gen, seed2)); case 3: - return extractMatPoly(gen_zonotope_exponential(dim_gen, m_gen)); - + return extractMatPoly(gen_zonotope_exponential(dim_gen, m_gen, seed2)); + } } else if (Vpoly_gen) { @@ -71,10 +76,10 @@ Rcpp::NumericMatrix poly_gen (int kind_gen, bool Vpoly_gen, bool Zono_gen, int d return extractMatPoly(gen_simplex(dim_gen, true)); case 4: - return extractMatPoly(random_vpoly(dim_gen, m_gen)); + return extractMatPoly(random_vpoly(dim_gen, m_gen, seed2)); case 5: - return extractMatPoly(random_vpoly_incube(dim_gen, m_gen)); + return extractMatPoly(random_vpoly_incube(dim_gen, m_gen, seed2)); } } else { @@ -96,7 +101,7 @@ Rcpp::NumericMatrix poly_gen (int kind_gen, bool Vpoly_gen, bool Zono_gen, int d return extractMatPoly(gen_skinny_cube(dim_gen)); case 6: - return extractMatPoly(random_hpoly(dim_gen, m_gen)); + return extractMatPoly(random_hpoly(dim_gen, m_gen, seed2)); } } diff --git a/R-proj/src/rotating.cpp b/R-proj/src/rotating.cpp index 4d5dbfdb2..002be1129 100644 --- a/R-proj/src/rotating.cpp +++ b/R-proj/src/rotating.cpp @@ -18,6 +18,7 @@ #include "hpolytope.h" #include "vpolytope.h" #include "zpolytope.h" +#include "vpolyintersectvpoly.h" #include "samplers.h" #include "rotating.h" #include "extractMatPoly.h" @@ -25,13 +26,14 @@ //' An internal Rccp function for the random rotation of a convex polytope //' //' @param P A convex polytope (H-, V-polytope or a zonotope). +//' @param T Optional. A rotation matrix. //' //' @section warning: //' Do not use this function. //' //' @return A matrix that describes the rotated polytope // [[Rcpp::export]] -Rcpp::NumericMatrix rotating (Rcpp::Reference P){ +Rcpp::NumericMatrix rotating (Rcpp::Reference P, Rcpp::Nullable T = R_NilValue){ typedef double NT; typedef Cartesian Kernel; @@ -40,13 +42,10 @@ Rcpp::NumericMatrix rotating (Rcpp::Reference P){ typedef HPolytope Hpolytope; typedef VPolytope Vpolytope; typedef Zonotope zonotope; + typedef IntersectionOfVpoly InterVP; typedef Eigen::Matrix VT; typedef Eigen::Matrix MT; - Hpolytope HP; - Vpolytope VP; - zonotope ZP; - MT TransorfMat; Rcpp::NumericMatrix Mat; unsigned int n = P.field("dimension"); @@ -57,7 +56,12 @@ Rcpp::NumericMatrix rotating (Rcpp::Reference P){ // Hpolytope Hpolytope HP; HP.init(n, Rcpp::as(P.field("A")), Rcpp::as(P.field("b"))); - TransorfMat = rotating < MT > (HP); + if (T.isNotNull()) { + TransorfMat = Rcpp::as(T); + HP.linear_transformIt(TransorfMat.inverse()); + } else { + TransorfMat = rotating < MT > (HP); + } Mat = extractMatPoly(HP); break; } @@ -65,7 +69,12 @@ Rcpp::NumericMatrix rotating (Rcpp::Reference P){ // Vpolytope Vpolytope VP; VP.init(n, Rcpp::as(P.field("V")), VT::Ones(Rcpp::as(P.field("V")).rows())); - TransorfMat = rotating< MT >(VP); + if (T.isNotNull()) { + TransorfMat = Rcpp::as(T); + VP.linear_transformIt(TransorfMat.inverse()); + } else { + TransorfMat = rotating < MT > (VP); + } Mat = extractMatPoly(VP); break; } @@ -73,12 +82,33 @@ Rcpp::NumericMatrix rotating (Rcpp::Reference P){ // Zonotope zonotope ZP; ZP.init(n, Rcpp::as(P.field("G")), VT::Ones(Rcpp::as(P.field("G")).rows())); - TransorfMat = rotating < MT > (ZP); + if (T.isNotNull()) { + TransorfMat = Rcpp::as(T); + ZP.linear_transformIt(TransorfMat.inverse()); + } else { + TransorfMat = rotating < MT > (ZP); + } Mat = extractMatPoly(ZP); break; } + case 4: { + Vpolytope VP1; + Vpolytope VP2; + InterVP VPcVP; + VP1.init(n, Rcpp::as(P.field("V1")), VT::Ones(Rcpp::as(P.field("V1")).rows())); + VP2.init(n, Rcpp::as(P.field("V2")), VT::Ones(Rcpp::as(P.field("V2")).rows())); + VPcVP.init(VP1, VP2); + if (T.isNotNull()) { + TransorfMat = Rcpp::as(T); + VPcVP.linear_transformIt(TransorfMat.inverse()); + } else { + TransorfMat = rotating < MT > (VPcVP); + } + } } + + TransorfMat.conservativeResize(n+1, n); TransorfMat.row(n) = VT::Ones(n); MT res(TransorfMat.rows(), Rcpp::as(Mat).rows()+n); diff --git a/R-proj/src/rounding.cpp b/R-proj/src/rounding.cpp index ae7f5f8e8..5e16d8481 100644 --- a/R-proj/src/rounding.cpp +++ b/R-proj/src/rounding.cpp @@ -27,19 +27,13 @@ //' Internal rcpp function for the rounding of a convex polytope //' //' @param P A convex polytope (H- or V-representation or zonotope). -//' @param random_walk Optional. A string that declares the random walk. -//' @param walk_length Optional. The number of the steps for the random walk. -//' @param parameters Optional. A list for the parameters of the methods: //' //' @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::Nullable random_walk = R_NilValue, - Rcpp::Nullable walk_length = R_NilValue, - Rcpp::Nullable parameters = R_NilValue){ +Rcpp::List rounding (Rcpp::Reference P){ typedef double NT; typedef Cartesian Kernel; @@ -61,7 +55,7 @@ Rcpp::List rounding (Rcpp::Reference P, NN=false, birk=false, verbose=false, - cdhr=true, rdhr = false, ball_walk = false, billiard = false; + cdhr=false, rdhr = false, ball_walk = false, billiard = false; NT delta = -1.0, diam = -1.0; unsigned int n = P.field("dimension"); @@ -72,28 +66,12 @@ Rcpp::List rounding (Rcpp::Reference P, Rcpp::NumericMatrix Mat; int type = P.field("type"); - if(!random_walk.isNotNull()) { - if (type == 1) { - cdhr = true; - } else { - billiard = true; - } - } else if(Rcpp::as(random_walk).compare(std::string("CDHR"))==0) { + if (type == 1) { + walkL = 10 + 10/n; cdhr = true; - } else if (Rcpp::as(random_walk).compare(std::string("RDHR"))==0) { - rdhr = true; - } else if (Rcpp::as(random_walk).compare(std::string("BW"))==0) { - if (Rcpp::as(parameters).containsElementNamed("BW_rad")) { - delta = Rcpp::as(Rcpp::as(parameters)["BW_rad"]); - } - ball_walk = true; - } else if (Rcpp::as(random_walk).compare(std::string("BiW")) == 0) { - billiard = true; - walkL = 1; - if (Rcpp::as(parameters).containsElementNamed("diameter")) - diam = Rcpp::as(Rcpp::as(parameters)["diameter"]); } else { - throw Rcpp::exception("Unknown walk type!"); + walkL = 5; + billiard = true; } switch (type) { @@ -101,20 +79,20 @@ 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); + //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(); - if (billiard && diam < 0.0) VP.comp_diam(diam, 0.0); + 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(); - if (billiard && diam < 0.0) ZP.comp_diam(diam, 0.0); + ZP.comp_diam(diam, 0.0); break; } case 4: { @@ -139,12 +117,6 @@ Rcpp::List rounding (Rcpp::Reference P, //default: throw Rcpp::exception("Wrong polytope input"); } - if(walk_length.isNotNull()) walkL = Rcpp::as(walk_length); - - if (ball_walk && delta < 0.0) { - delta = 4.0 * InnerBall.second / std::sqrt(NT(n)); - } - unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); // the random engine with this seed RNGType rng(seed); diff --git a/R-proj/src/sample_points.cpp b/R-proj/src/sample_points.cpp index 5ed69aae7..d77c16f2e 100644 --- a/R-proj/src/sample_points.cpp +++ b/R-proj/src/sample_points.cpp @@ -2,10 +2,10 @@ // VolEsti (volume computation and sampling library) -// Copyright (c) 20012-2018 Vissarion Fisikopoulos -// Copyright (c) 2018 Apostolos Chalkis +// Copyright (c) 20012-2020 Vissarion Fisikopoulos +// Copyright (c) 2018-2020 Apostolos Chalkis -//Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 program. +//Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 and 2019 program. #include #include @@ -31,22 +31,27 @@ //' Sample n points with uniform or multidimensional spherical gaussian -centered in an internal point- target distribution. //' The \eqn{d}-dimensional unit simplex is the set of points \eqn{\vec{x}\in \R^d}, s.t.: \eqn{\sum_i x_i\leq 1}, \eqn{x_i\geq 0}. The \eqn{d}-dimensional canonical simplex is the set of points \eqn{\vec{x}\in \R^d}, s.t.: \eqn{\sum_i x_i = 1}, \eqn{x_i\geq 0}. //' -//' @param P A convex polytope. It is an object from class (a) Hpolytope or (b) Vpolytope or (c) Zonotope. -//' @param n The number of points that the function is going to sample from the convex polytope. The default value is \eqn{100}. -//' @param distribution Optional. A string that declares the target distribution: a) \code{'uniform'} for the uniform distribution or b) \code{'gaussian'} for the multidimensional spherical distribution. The default target distribution is uniform. -//' @param random_walk Optional. 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 d) \code{'BiW'} for Billiard walk. The default walk is \code{'BiW'} for the uniform distribution or \code{'CDHR'} for the Normal distribution. -//' @param walk_length Optional. The number of the steps for the random walk. The default value is \eqn{1} for \code{'BiW'} and \eqn{\lfloor 10 + d/10\rfloor} otherwise, where \eqn{d} is the dimension that the polytope lies. -//' @param exact A boolean parameter. It should be used for the uniform sampling from the boundary or the interior of a hypersphere centered at the origin or from the unit or the canonical or an arbitrary simplex. The arbitrary simplex has to be given as a V-polytope. For the rest well known convex bodies the dimension has to be declared and the type of body as well as the radius of the hypersphere. -//' @param body A string that declares the type of the body for the exact sampling: a) \code{'unit simplex'} for the unit simplex, b) \code{'canonical simplex'} for the canonical simplex, c) \code{'hypersphere'} for the boundary of a hypersphere centered at the origin, d) \code{'ball'} for the interior of a hypersphere centered at the origin. -//' @param parameters Optional. A list for the parameters of the methods: +//' @param P A convex polytope. It is an object from class (a) Hpolytope or (b) Vpolytope or (c) Zonotope or (d) an intersection of two V-polytopes. +//' @param n The number of points that the function is going to sample from the convex polytope. +//' @param random_walk Optional. A list that declares the random walk and some related parameters as follows: //' \itemize{ +//' \item{\code{walk} }{ A string to declare the random walk: i) \code{'CDHR'} for Coordinate Directions Hit-and-Run, ii) \code{'RDHR'} for Random Directions Hit-and-Run, iii) \code{'BaW'} for Ball Walk, iv) \code{'BiW'} for Billiard walk, v) \code{'BCDHR'} boundary sampling by keeping the extreme points of CDHR or vi) \code{'BRDHR'} boundary sampling by keeping the extreme points of RDHR. The default walk is \code{'BiW'} for the uniform distribution or \code{'CDHR'} for the Gaussian distribution.} +//' \item{\code{walk_length} }{ The number of the steps for the random walk. The default value is \eqn{5} for \code{'BiW'} and \eqn{\lfloor 10 + d/10\rfloor} otherwise.} +//' \item{\code{BaW_rad} }{ The radius for the ball walk.} +//' \item{\code{L} }{The maximum length of the billiard trajectory.} +//' } +//' @param distribution Optional. A list that declares the target density and some related parameters as follows: +//' \itemize{ +//' \item{\code{density}}{A string: (a) \code{'uniform'} for the uniform distribution or b) \code{'gaussian'} for the multidimensional spherical distribution. The default target distribution is uniform.} //' \item{\code{variance} }{ The variance of the multidimensional spherical gaussian. The default value is 1.} +//' \item{\code{InnerPoint} }{ A \eqn{d}-dimensional numerical vector that defines a starting point in the interior of the polytope for the random walk and the mode of the Gaussian distribution. The default choice is the center of the Chebychev ball.} +//' } +//' @param known_body A list to request exact uniform sampling from special well known convex bodies through the following input parameters: +//' \itemize{ +//' \item{\code{body} }{ A string that declares the type of the body for the exact sampling: a) \code{'unit simplex'} for the unit simplex, b) \code{'canonical simplex'} for the canonical simplex, c) \code{'hypersphere'} for the boundary of a hypersphere centered at the origin, d) \code{'ball'} for the interior of a hypersphere centered at the origin.} //' \item{\code{dimension} }{ An integer that declares the dimension when exact sampling is enabled for a simplex or a hypersphere.} //' \item{\code{radius} }{ The radius of the \eqn{d}-dimensional hypersphere. The default value is \eqn{1}.} -//' \item{\code{BW_rad} }{ The radius for the ball walk.} -//' \item{\code{L} }{The maximum length of the billiard trajectory.} //' } -//' @param InnerPoint A \eqn{d}-dimensional numerical vector that defines a point in the interior of polytope P. //' //' @references \cite{R.Y. Rubinstein and B. Melamed, //' \dQuote{Modern simulation and modeling} \emph{ Wiley Series in Probability and Statistics,} 1998.} @@ -59,32 +64,27 @@ //' @examples //' # uniform distribution from the 3d unit cube in V-representation using ball walk //' P = gen_cube(3, 'V') -//' points = sample_points(P, random_walk = "BaW", walk_length = 5) +//' points = sample_points(P, n = 100, random_walk = list("walk" = "BaW", "walk_length" = 5)) //' //' # gaussian distribution from the 2d unit simplex in H-representation with variance = 2 //' A = matrix(c(-1,0,0,-1,1,1), ncol=2, nrow=3, byrow=TRUE) //' b = c(0,0,1) //' P = Hpolytope$new(A,b) -//' points = sample_points(P, distribution = "gaussian", parameters = list("variance" = 2)) +//' points = sample_points(P, n = 100, distribution = list("density" = "gaussian", "variance" = 2)) //' -//' # uniform points from the boundary of a 10-dimensional hypersphere -//' points = sample_points(exact = TRUE, body = "hypersphere", parameters = list("dimension" = 10)) +//' # uniform points from the boundary of a 2-dimensional random H-polytope +//' P = gen_rand_hpoly(2,20) +//' points = sample_points(P, n = 5000, random_walk = list("walk" = "BRDHR")) //' -//' # 10000 uniform points from a 2-d arbitrary simplex -//' V = matrix(c(2,3,-1,7,0,0),ncol = 2, nrow = 3, byrow = TRUE) -//' P = Vpolytope$new(V) -//' points = sample_points(P, n = 10000, exact = TRUE) +//' # 100 uniform points from the 2-d unit ball +//' points = sample_points(n = 100, known_body = list("body" = "ball", "dimension" = 2)) //' @export // [[Rcpp::export]] 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 random_walk = R_NilValue, - Rcpp::Nullable walk_length = R_NilValue, - Rcpp::Nullable exact = R_NilValue, - Rcpp::Nullable body = R_NilValue, - Rcpp::Nullable parameters = R_NilValue, - Rcpp::Nullable InnerPoint = R_NilValue){ + Rcpp::Nullable known_body = R_NilValue){ typedef double NT; typedef Cartesian Kernel; @@ -104,62 +104,58 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P = R_NilValue int type, dim, numpoints; NT radius = 1.0, delta = -1.0, diam = -1.0; - bool set_mean_point = false, cdhr = false, rdhr = false, ball_walk = false, gaussian = false, billiard = false; + bool set_mean_point = false, cdhr = false, rdhr = false, ball_walk = false, gaussian = false, + billiard = false, boundary = false; std::list randPoints; std::pair InnerBall; + Point shift(dim); numpoints = (!n.isNotNull()) ? 100 : Rcpp::as(n); - if (exact.isNotNull()) { + if (!n.isNotNull()) { + throw Rcpp::exception("The number of samples is not declared!"); + } else { + numpoints = Rcpp::as(n); + if (numpoints <= 0) throw Rcpp::exception("The number of samples has to be a positice integer!"); + } + + if (known_body.isNotNull()) { if (P.isNotNull()) { - type = Rcpp::as(P).field("type"); - dim = Rcpp::as(P).field("dimension"); - if (Rcpp::as(exact) && type==2) { - if (Rcpp::as(Rcpp::as(P).field("V")).rows() == Rcpp::as(Rcpp::as(P).field("V")).cols()+1) { - Vpolytope VP; - VP.init(dim, - Rcpp::as(Rcpp::as(P).field("V")), - VT::Ones(Rcpp::as(Rcpp::as(P).field("V")).rows())); - Sam_arb_simplex(VP, numpoints, randPoints); - } else { - throw Rcpp::exception("Not a simplex!"); - } - } else if (Rcpp::as(exact) && type!=2) { - throw Rcpp::exception("Not a simplex in V-representation!"); - } + throw Rcpp::exception("No input Polytope is necessary when a known body is declared!"); } else { - if (!body.isNotNull()) { + if (!Rcpp::as(known_body).containsElementNamed("dimension")) { - throw Rcpp::exception("Wrong input!"); + throw Rcpp::exception("Dimension has to be given as input!"); - } else if (!Rcpp::as(parameters).containsElementNamed("dimension")) { + } + dim = Rcpp::as(Rcpp::as(known_body)["dimension"]); + if (Rcpp::as(known_body).containsElementNamed("radius")) { - throw Rcpp::exception("Wrong input!"); + radius = Rcpp::as(Rcpp::as(known_body)["radius"]); } - dim = Rcpp::as(Rcpp::as(parameters)["dimension"]); - if (Rcpp::as(parameters).containsElementNamed("radius")) { + if (!Rcpp::as(known_body).containsElementNamed("body")) { - radius = Rcpp::as(Rcpp::as(parameters)["radius"]); + throw Rcpp::exception("The kind of body has to be given as input!"); } - if (Rcpp::as(body).compare(std::string("hypersphere"))==0) { + if (Rcpp::as(Rcpp::as(known_body)["body"]).compare(std::string("hypersphere"))==0) { for (unsigned int k = 0; k < numpoints; ++k) { randPoints.push_back(get_point_on_Dsphere(dim, radius)); } - } else if (Rcpp::as(body).compare(std::string("ball"))==0) { + } else if (Rcpp::as(Rcpp::as(known_body)["body"]).compare(std::string("ball"))==0) { for (unsigned int k = 0; k < numpoints; ++k) { randPoints.push_back(get_point_in_Dsphere(dim, radius)); } - } else if (Rcpp::as(body).compare(std::string("unit simplex"))==0) { + } else if (Rcpp::as(Rcpp::as(known_body)["body"]).compare(std::string("unit simplex"))==0) { Sam_Unit(dim, numpoints, randPoints); - } else if (Rcpp::as(body).compare(std::string("canonical simplex"))==0) { + } else if (Rcpp::as(Rcpp::as(known_body)["body"]).compare(std::string("canonical simplex"))==0) { Sam_Canon_Unit(dim, numpoints, randPoints); @@ -175,56 +171,85 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P = R_NilValue dim = Rcpp::as(P).field("dimension"); unsigned int walkL = 10 + dim / 10; - if (!distribution.isNotNull() || Rcpp::as(distribution).compare(std::string("uniform")) == 0) { - walkL = 3; + //std::cout<< Rcpp::as(distribution).containsElementNamed("density") <(Rcpp::as(distribution)["density"]).compare(std::string("uniform"))<(distribution).containsElementNamed("density")) { billiard = true; - } else if (Rcpp::as(distribution).compare(std::string("gaussian")) == 0) { + } else if ( + Rcpp::as(Rcpp::as(distribution)["density"]).compare(std::string("uniform")) == + 0) { + billiard = true; + } else if ( + Rcpp::as(Rcpp::as(distribution)["density"]).compare(std::string("gaussian")) == + 0) { gaussian = true; - cdhr = true; } else { throw Rcpp::exception("Wrong distribution!"); } Point MeanPoint; - if (InnerPoint.isNotNull()) { - if (Rcpp::as(InnerPoint).size() != dim) { - Rcpp::warning("Internal Point has to lie in the same dimension as the polytope P"); + if (Rcpp::as(distribution).containsElementNamed("inner_point")) { + if (Rcpp::as(Rcpp::as(distribution)["inner_point"]).size() != dim) { + throw Rcpp::exception("Internal Point has to lie in the same dimension as the polytope P"); } else { set_mean_point = true; - MeanPoint = Point(dim, Rcpp::as < std::vector < NT > > (InnerPoint).begin(), - Rcpp::as < std::vector < NT > > (InnerPoint).end()); + VT temp = Rcpp::as(Rcpp::as(distribution)["inner_point"]); + MeanPoint = Point(dim, std::vector(&temp[0], temp.data() + temp.cols() * temp.rows())); + //MeanPoint = Point(dim, Rcpp::as < std::vector < NT > > (Rcpp::as(distribution)["inner_point"]).begin(), + // Rcpp::as < std::vector < NT > > (Rcpp::as(distribution)["inner_point"]).end()); } } - if (walk_length.isNotNull()) walkL = Rcpp::as(walk_length); NT a = 0.5; - - if (Rcpp::as(parameters).containsElementNamed("variance")) - a = 1.0 / (2.0 * Rcpp::as(Rcpp::as(parameters)["variance"])); + if (Rcpp::as(distribution).containsElementNamed("variance")) { + a = 1.0 / (2.0 * Rcpp::as(Rcpp::as(distribution)["variance"])); + if (!gaussian) { + Rcpp::warning("The variance can be set only for Gaussian sampling!"); + } else if (a <= 0.0) { + throw Rcpp::exception("The variance has to be positive!"); + } + } if (!random_walk.isNotNull()) { if (gaussian) { cdhr = true; } else { billiard = true; + walkL = 5; } - } else if (Rcpp::as(random_walk).compare(std::string("CDHR")) == 0) { + } else if (Rcpp::as(Rcpp::as(random_walk)["walk"]).compare(std::string("CDHR")) == 0) { cdhr = true; - } else if (Rcpp::as(random_walk).compare(std::string("RDHR")) == 0) { + } else if (Rcpp::as(Rcpp::as(random_walk)["walk"]).compare(std::string("RDHR")) == 0) { rdhr = true; - } else if (Rcpp::as(random_walk).compare(std::string("BaW")) == 0) { - if (Rcpp::as(parameters).containsElementNamed("BW_rad")) - delta = Rcpp::as(Rcpp::as(parameters)["BW_rad"]); + } 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"]); ball_walk = true; - } else if (Rcpp::as(random_walk).compare(std::string("BiW")) == 0) { + } else if (Rcpp::as(Rcpp::as(random_walk)["walk"]).compare(std::string("BiW")) == 0) { if (gaussian) throw Rcpp::exception("Billiard walk can be used only for uniform sampling!"); billiard = true; - if (Rcpp::as(parameters).containsElementNamed("diameter")) - diam = Rcpp::as(Rcpp::as(parameters)["diameter"]); - } else { + walkL = 5; + if (Rcpp::as(random_walk).containsElementNamed("L")) + diam = Rcpp::as(Rcpp::as(random_walk)["L"]); + } 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!"); + rdhr = true; + boundary = true; + }else if (Rcpp::as(Rcpp::as(random_walk)["walk"]).compare(std::string("BCDHR")) == 0) { + if (gaussian) throw Rcpp::exception("Gaussian sampling from the boundary is not supported!"); + cdhr = true; + boundary = true; + }else { throw Rcpp::exception("Unknown walk type!"); } + if (Rcpp::as(random_walk).containsElementNamed("walk_length")) { + walkL = Rcpp::as(Rcpp::as(random_walk)["walk_length"]); + if (walkL <= 0) { + throw Rcpp::exception("The walk length has to be a positive integer!"); + } + } + bool rand_only=false, NN=false, birk=false, @@ -246,7 +271,15 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P = R_NilValue InnerBall = HP.ComputeInnerBall(); if (!set_mean_point) MeanPoint = InnerBall.first; } + if (HP.is_in(MeanPoint) == 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) { + shift = MeanPoint; + HP.shift(Eigen::Map(&MeanPoint.get_coeffs()[0], MeanPoint.dimension())); + MeanPoint = Point(dim); + } break; } case 2: { @@ -258,7 +291,14 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P = R_NilValue InnerBall = VP.ComputeInnerBall(); if (!set_mean_point) MeanPoint = InnerBall.first; } + if (VP.is_in(MeanPoint) == 0) + throw Rcpp::exception("The given point is not in the interior of the polytope!"); if (billiard && diam < 0.0) VP.comp_diam(diam, 0.0); + if (gaussian) { + shift = MeanPoint; + VP.shift(Eigen::Map(&MeanPoint.get_coeffs()[0], MeanPoint.dimension())); + MeanPoint = Point(dim); + } break; } case 3: { @@ -270,7 +310,14 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P = R_NilValue InnerBall = ZP.ComputeInnerBall(); if (!set_mean_point) MeanPoint = InnerBall.first; } + if (ZP.is_in(MeanPoint) == 0) + throw Rcpp::exception("The given point is not in the interior of the polytope!"); if (billiard && diam < 0.0) ZP.comp_diam(diam, 0.0); + if (gaussian) { + shift = MeanPoint; + ZP.shift(Eigen::Map(&MeanPoint.get_coeffs()[0], MeanPoint.dimension())); + MeanPoint = Point(dim); + } break; } case 4: { @@ -286,9 +333,16 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P = R_NilValue if (!VPcVP.is_feasible()) throw Rcpp::exception("Empty set!"); InnerBall = VPcVP.ComputeInnerBall(); if (!set_mean_point) MeanPoint = InnerBall.first; + if (VPcVP.is_in(MeanPoint) == 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) { + shift = MeanPoint; + VPcVP.shift(Eigen::Map(&MeanPoint.get_coeffs()[0], MeanPoint.dimension())); + MeanPoint = Point(dim); + } break; } } @@ -309,37 +363,43 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P = R_NilValue switch (type) { case 1: { sampling_only(randPoints, HP, walkL, numpoints, gaussian, - a, MeanPoint, var1, var2); + a, boundary, MeanPoint, var1, var2); break; } case 2: { sampling_only(randPoints, VP, walkL, numpoints, gaussian, - a, MeanPoint, var1, var2); + a, boundary, MeanPoint, var1, var2); break; } case 3: { sampling_only(randPoints, ZP, walkL, numpoints, gaussian, - a, MeanPoint, var1, var2); + a, boundary, MeanPoint, var1, var2); break; } case 4: { sampling_only(randPoints, VPcVP, walkL, numpoints, gaussian, - a, MeanPoint, var1, var2); + a, boundary, MeanPoint, var1, var2); break; } } } else { - throw Rcpp::exception("Wrong input!"); + throw Rcpp::exception("Wrong inputs!"); } + if (numpoints % 2 == 1 && boundary) numpoints--; MT RetMat(dim, numpoints); unsigned int jj = 0; - for (typename std::list::iterator rpit = randPoints.begin(); rpit!=randPoints.end(); rpit++, jj++) - RetMat.col(jj) = Eigen::Map(&(*rpit).get_coeffs()[0], (*rpit).dimension()); + for (typename std::list::iterator rpit = randPoints.begin(); rpit!=randPoints.end(); rpit++, jj++) { + if (gaussian) { + RetMat.col(jj) = Eigen::Map(&(*rpit).get_coeffs()[0], (*rpit).dimension()) + Eigen::Map(&shift.get_coeffs()[0], shift.dimension()); + } else { + RetMat.col(jj) = Eigen::Map(&(*rpit).get_coeffs()[0], (*rpit).dimension()); + } + } return Rcpp::wrap(RetMat); } diff --git a/R-proj/src/volume.cpp b/R-proj/src/volume.cpp index 2414f0e6e..814050658 100644 --- a/R-proj/src/volume.cpp +++ b/R-proj/src/volume.cpp @@ -2,10 +2,10 @@ // VolEsti (volume computation and sampling library) -// Copyright (c) 20012-2018 Vissarion Fisikopoulos -// Copyright (c) 2018 Apostolos Chalkis +// Copyright (c) 20012-2020 Vissarion Fisikopoulos +// Copyright (c) 2018-2020 Apostolos Chalkis -//Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 program. +//Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 and 2019 program. #include #include @@ -20,7 +20,7 @@ template double generic_volume(Polytope& P, unsigned int walk_step, double e, - Rcpp::Nullable InnerBall, bool CG, bool CB, bool hpoly, unsigned int win_len, + 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) @@ -45,16 +45,10 @@ double generic_volume(Polytope& P, unsigned int walk_step, double e, std::pair InnerB; - if(InnerBall.isNotNull()) { //if it is given as an input - // store internal point hat is given as input - Rcpp::NumericVector InnerVec = Rcpp::as(InnerBall); - std::vector temp_p; - for (unsigned int j=0; j 0.0) { //if it is given as an input + + InnerB = InnerBall; + } else if(type == 2 && CB) { if (rounding) { @@ -125,30 +119,30 @@ double generic_volume(Polytope& P, unsigned int walk_step, double e, //' For the volume approximation can be used two algorithms. Either SequenceOfBalls or CoolingGaussian. A H-polytope with \eqn{m} facets is described by a \eqn{m\times d} matrix \eqn{A} and a \eqn{m}-dimensional vector \eqn{b}, s.t.: \eqn{Ax\leq b}. A V-polytope is defined as the convex hull of \eqn{m} \eqn{d}-dimensional points which correspond to the vertices of P. A zonotope is desrcibed by the Minkowski sum of \eqn{m} \eqn{d}-dimensional segments. //' //' @param P A convex polytope. It is an object from class (a) Hpolytope or (b) Vpolytope or (c) Zonotope. -//' @param walk_length Optional. The number of the steps for the random walk. The default value is \eqn{\lfloor 10 + d/10\rfloor} for SequenceOfBalls and \eqn{1} otherwise. -//' @param error Optional. Declare the upper bound for the approximation error. The default value is \eqn{1} for SequenceOfBalls and \eqn{0.1} otherwise. -//' @param inner_ball Optional. A \eqn{d+1} 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. -//' @param algo Optional. A string that declares which algorithm to use: a) \code{'SoB'} for SequenceOfBalls or b) \code{'CG'} for CoolingGaussian or c) \code{'CB'} for cooling bodies. -//' @param random_walk Optional. 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. -//' @param rounding Optional. A boolean parameter for rounding. The default value is \code{TRUE} for V-polytopes and \code{FALSE} otherwise. -//' @param parameters Optional. A list for the parameters of the algorithms: +//' @param algo Optional. A list that declares which algorithm, random walk and values of parameters to use, as follows: //' \itemize{ -//' \item{\code{Window} }{ 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{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 stopping criterion.} -//' \item{\code{L} }{The maximum length of the billiard trajectory.} +//' \item{\code{algorithm} }{ A string to set the algorithm to use: a) \code{'SoB'} for SequenceOfBalls or b) \code{'CG'} for CoolingGaussian or c) \code{'CB'} for cooling bodies. The defalut algorithm for H-polytopes is \code{'CB'} when \eqn{d\leq 200} and \code{'CG'} when \eqn{d>200}. For the other representations the default algorithm is \code{'CB'}.} +//' \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{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. //' //' @references \cite{I.Z.Emiris and V. Fisikopoulos, //' \dQuote{Practical polytope volume approximation,} \emph{ACM Trans. Math. Soft.,} 2014.}, @@ -165,20 +159,16 @@ double generic_volume(Polytope& P, unsigned int walk_step, double e, //' //' # calling CG algorithm for a V-polytope (3d simplex) //' P = gen_simplex(2,'V') -//' vol = volume(P, algo = "CG") +//' vol = volume(P, algo = list("algorithm" = "CG")) //' //' # calling CG algorithm for a 2-dimensional zonotope defined as the Minkowski sum of 4 segments //' Z = gen_rand_zonotope(2, 4) -//' vol = volume(Z, random_walk = "RDHR", walk_length = 5) +//' vol = volume(Z, algo = list("random_walk" = "RDHR", "walk_length" = 5)) //' @export // [[Rcpp::export]] -double volume (Rcpp::Reference P, Rcpp::Nullable walk_length = R_NilValue, - Rcpp::Nullable error = R_NilValue, - Rcpp::Nullable inner_ball = R_NilValue, - Rcpp::Nullable algo = R_NilValue, - Rcpp::Nullable random_walk = R_NilValue, - Rcpp::Nullable rounding = R_NilValue, - Rcpp::Nullable parameters = R_NilValue) { +double volume (Rcpp::Reference P, + Rcpp::Nullable algo = R_NilValue, + Rcpp::Nullable rounding = R_NilValue) { typedef double NT; typedef Cartesian Kernel; @@ -193,47 +183,54 @@ double volume (Rcpp::Reference P, Rcpp::Nullable walk_length = R_ unsigned int n = P.field("dimension"), walkL; int type = P.field("type"); - bool CG = false, CB = false, cdhr = false, rdhr = false, ball_walk = false, round = false, win2 = false, hpoly = false, - billiard = false; - unsigned int win_len = 4*n*n+500, N = 500 * 2 + n * n / 2, NN = 120 + (n*n)/10, nu = 10; + 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; 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; - if(!algo.isNotNull()){ - + if (!Rcpp::as(algo).containsElementNamed("algorithm")) { if (type == 2 || type == 3) { CB = true; - } else if (n<=200) { + } else if (n <= 200) { CB = true; } else { CG = true; } - e = (!error.isNotNull()) ? 0.1 : Rcpp::as(error); - walkL = (!walk_length.isNotNull()) ? 1 : Rcpp::as(walk_length); + walkL = (!Rcpp::as(algo).containsElementNamed("walk_length")) ? 1 : Rcpp::as( + Rcpp::as(algo)["walk_length"]); + e = (!Rcpp::as(algo).containsElementNamed("error")) ? 0.1 : Rcpp::as( + Rcpp::as(algo)["error"]); + } else if (Rcpp::as(Rcpp::as(algo)["algorithm"]).compare(std::string("SOB")) == 0) { - }else if (Rcpp::as(algo).compare(std::string("SOB"))==0){ + walkL = (!Rcpp::as(algo).containsElementNamed("walk_length")) ? 10 + n / 10 : Rcpp::as( + Rcpp::as(algo)["walk_length"]); + e = (!Rcpp::as(algo).containsElementNamed("error")) ? 1.0 : Rcpp::as( + Rcpp::as(algo)["error"]); - walkL = (!walk_length.isNotNull()) ? 10 + n / 10 : Rcpp::as(walk_length); - e = (!error.isNotNull()) ? 1.0 : Rcpp::as(error); - - } else if (Rcpp::as(algo).compare(std::string("CG"))==0) { + } else if (Rcpp::as(Rcpp::as(algo)["algorithm"]).compare(std::string("CG")) == 0) { CG = true; - e = (!error.isNotNull()) ? 0.1 : Rcpp::as(error); - walkL = (!walk_length.isNotNull()) ? 1 : Rcpp::as(walk_length); + walkL = (!Rcpp::as(algo).containsElementNamed("walk_length")) ? 1 : Rcpp::as( + Rcpp::as(algo)["walk_length"]); + e = (!Rcpp::as(algo).containsElementNamed("error")) ? 0.1 : Rcpp::as( + Rcpp::as(algo)["error"]); - } else if (Rcpp::as(algo).compare(std::string("CB")) == 0) { + } else if (Rcpp::as(Rcpp::as(algo)["algorithm"]).compare(std::string("CB")) == 0) { CB = true; - e = (!error.isNotNull()) ? 0.1 : Rcpp::as(error); - walkL = (!walk_length.isNotNull()) ? 1 : Rcpp::as(walk_length); + walkL = (!Rcpp::as(algo).containsElementNamed("walk_length")) ? 1 : Rcpp::as( + Rcpp::as(algo)["walk_length"]); + e = (!Rcpp::as(algo).containsElementNamed("error")) ? 0.1 : Rcpp::as( + Rcpp::as(algo)["error"]); } else { throw Rcpp::exception("Unknown method!"); } - if (!random_walk.isNotNull()) { + + if (!Rcpp::as(algo).containsElementNamed("random_walk")) { if ( type == 1 ){ cdhr = true; if (CB) win_len = 3*n*n+400; @@ -246,14 +243,14 @@ double volume (Rcpp::Reference P, Rcpp::Nullable walk_length = R_ rdhr = true; } } - }else if (Rcpp::as(random_walk).compare(std::string("CDHR")) == 0) { + }else if (Rcpp::as(Rcpp::as(algo)["random_walk"]).compare(std::string("CDHR")) == 0) { cdhr = true; if (CB) win_len = 3*n*n+400; - } else if (Rcpp::as(random_walk).compare(std::string("RDHR")) == 0) { + } else if (Rcpp::as(Rcpp::as(algo)["random_walk"]).compare(std::string("RDHR")) == 0) { rdhr = true; - } else if (Rcpp::as(random_walk).compare(std::string("BaW"))==0) { + } else if (Rcpp::as(Rcpp::as(algo)["random_walk"]).compare(std::string("BaW")) == 0) { ball_walk = true; - } else if (Rcpp::as(random_walk).compare(std::string("BiW"))==0) { + } else if (Rcpp::as(Rcpp::as(algo)["random_walk"]).compare(std::string("BiW")) == 0) { if (CG) { if (type !=1){ Rcpp::Rcout << "Billiard walk is not supported for CG algorithm. RDHR is used."< walk_length = R_ round = (!rounding.isNotNull()) ? false : Rcpp::as(rounding); } - if(parameters.isNotNull()) { + if (Rcpp::as(algo).containsElementNamed("BW_rad")) { + delta = Rcpp::as(Rcpp::as(algo)["BW_rad"]); + } + if (Rcpp::as(algo).containsElementNamed("C")) { + C = Rcpp::as(Rcpp::as(algo)["C"]); + N = 500 * ((int) C) + n * n / 2; + cg_params++; + } + if (Rcpp::as(algo).containsElementNamed("M")) { + N = Rcpp::as(Rcpp::as(algo)["M"]); + cg_params++; + } + if (Rcpp::as(algo).containsElementNamed("len_win")) { + win_len = Rcpp::as(Rcpp::as(algo)["len_win"]); + if (!CB && !CG) Rf_warning("flag 'len_win' can be used only for CG or CB algorithms."); + } + if (Rcpp::as(algo).containsElementNamed("frac")) { + frac = Rcpp::as(Rcpp::as(algo)["frac"]); + cg_params++; + } + if (Rcpp::as(algo).containsElementNamed("ratio")) { + ratio = Rcpp::as(Rcpp::as(algo)["ratio"]); + cg_params++; + } + if (Rcpp::as(algo).containsElementNamed("hpoly")) { + hpoly = Rcpp::as(Rcpp::as(algo)["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(algo).containsElementNamed("lb")) { + lb = Rcpp::as(Rcpp::as(algo)["lb"]); + cb_params++; + } + if (Rcpp::as(algo).containsElementNamed("ub")) { + ub = Rcpp::as(Rcpp::as(algo)["ub"]); + cb_params++; + } + if (Rcpp::as(algo).containsElementNamed("nu")) { + nu = Rcpp::as(Rcpp::as(algo)["nu"]); + cb_params++; + } + if (Rcpp::as(algo).containsElementNamed("N")) { + NN = Rcpp::as(Rcpp::as(algo)["N"]); + cb_params++; + } + if (Rcpp::as(algo).containsElementNamed("minmaxW")) { + win2 = Rcpp::as(Rcpp::as(algo)["minmaxW"]); + cb_params++; + } + if (Rcpp::as(algo).containsElementNamed("prob")) { + p = Rcpp::as(Rcpp::as(algo)["prob"]); + cb_params++; + } + if (Rcpp::as(algo).containsElementNamed("alpha")) { + alpha = Rcpp::as(Rcpp::as(algo)["alpha"]); + cb_params++; + } + if (Rcpp::as(algo).containsElementNamed("L")) { + diam = Rcpp::as(Rcpp::as(algo)["L"]); + cb_params++; + } + + 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."); + } - if (Rcpp::as(parameters).containsElementNamed("BW_rad")) { - delta = Rcpp::as(Rcpp::as(parameters)["BW_rad"]); - } - if (Rcpp::as(parameters).containsElementNamed("C")) { - C = Rcpp::as(Rcpp::as(parameters)["C"]); - N = 500 * ((int) C) + n * n / 2; - } - if (Rcpp::as(parameters).containsElementNamed("M")) { - N = Rcpp::as(Rcpp::as(parameters)["M"]); - } - if (Rcpp::as(parameters).containsElementNamed("Window")) { - win_len = Rcpp::as(Rcpp::as(parameters)["Window"]); - } - if (Rcpp::as(parameters).containsElementNamed("frac")) { - frac = Rcpp::as(Rcpp::as(parameters)["frac"]); - } - if (Rcpp::as(parameters).containsElementNamed("ratio")) { - ratio = Rcpp::as(Rcpp::as(parameters)["ratio"]); - } - if (Rcpp::as(parameters).containsElementNamed("hpoly")) { - hpoly = Rcpp::as(Rcpp::as(parameters)["hpoly"]); - 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(parameters).containsElementNamed("lb")) { - lb = Rcpp::as(Rcpp::as(parameters)["lb"]); - } - if (Rcpp::as(parameters).containsElementNamed("ub")) { - ub = Rcpp::as(Rcpp::as(parameters)["ub"]); - } - if (Rcpp::as(parameters).containsElementNamed("nu")) { - nu = Rcpp::as(Rcpp::as(parameters)["nu"]); - } - if (Rcpp::as(parameters).containsElementNamed("N")) { - NN = Rcpp::as(Rcpp::as(parameters)["N"]); - } - if (Rcpp::as(parameters).containsElementNamed("minmaxW")) { - win2 = Rcpp::as(Rcpp::as(parameters)["minmaxW"]); - } - if (Rcpp::as(parameters).containsElementNamed("prob")) { - p = Rcpp::as(Rcpp::as(parameters)["prob"]); - } - if (Rcpp::as(parameters).containsElementNamed("alpha")) { - alpha = Rcpp::as(Rcpp::as(parameters)["alpha"]); - } - if (Rcpp::as(parameters).containsElementNamed("diameter")) { - diam = Rcpp::as(Rcpp::as(parameters)["diameter"]); + std::pair inner_ball; + inner_ball.second = -1.0; + if (Rcpp::as(algo).containsElementNamed("inner_ball")) { + if (Rcpp::as(Rcpp::as(algo)["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(algo)["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."); } } @@ -346,6 +372,9 @@ double volume (Rcpp::Reference P, Rcpp::Nullable walk_length = R_ // 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); } @@ -353,6 +382,9 @@ double volume (Rcpp::Reference P, Rcpp::Nullable walk_length = R_ // 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); } @@ -360,6 +392,9 @@ double volume (Rcpp::Reference P, Rcpp::Nullable walk_length = R_ // 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!"); + } 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); } @@ -371,8 +406,10 @@ double volume (Rcpp::Reference P, Rcpp::Nullable walk_length = R_ VP1.init(n, Rcpp::as(P.field("V1")), VT::Ones(Rcpp::as(P.field("V1")).rows())); VP2.init(n, Rcpp::as(P.field("V2")), VT::Ones(Rcpp::as(P.field("V2")).rows())); VPcVP.init(VP1, VP2); - Rcpp::NumericVector InnerVec(n + 1); 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); } diff --git a/R-proj/src/zonotope_approximation.cpp b/R-proj/src/zonotope_approximation.cpp index b2961ba92..a34b9dfb5 100644 --- a/R-proj/src/zonotope_approximation.cpp +++ b/R-proj/src/zonotope_approximation.cpp @@ -19,25 +19,16 @@ //' //' @param Z A zonotope. //' @param fit_ratio Optional. A boolean parameter to request the computation of the ratio of fitness. -//' @param walk_length Optional. The number of the steps for the random walk. The default value is \eqn{\lfloor 10 + d/10\rfloor} for SequenceOfBalls and \eqn{1} for CoolingGaussian. -//' @param error Optional. Declare the upper bound for the approximation error. The default value is \eqn{1} for SequenceOfBalls and \eqn{0.1} for CoolingGaussian. -//' @param inner_ball Optional. A \eqn{d+1} 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. -//' @param random_walk Optional. 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 or c) \code{'BW'} for Ball Walk. The default walk is \code{'CDHR'}. -//' @param rounding Optional. A boolean parameter for rounding. The default value is \code{FALSE}. -//' @param parameters Optional. A list for the parameters of the volume algorithm +//' @param algo_parameters Optional. A list that declares the values of the parameters of CB algorithm. //' //' @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. // [[Rcpp::export]] -Rcpp::List zono_approx (Rcpp::Reference Z, Rcpp::Nullable fit_ratio = R_NilValue, - Rcpp::Nullable walk_length = R_NilValue, - Rcpp::Nullable error = R_NilValue, - Rcpp::Nullable inner_ball = R_NilValue, - Rcpp::Nullable random_walk = R_NilValue, - Rcpp::Nullable rounding = R_NilValue, - Rcpp::Nullable parameters = R_NilValue) { +Rcpp::List zono_approx (Rcpp::Reference Z, + Rcpp::Nullable fit_ratio = R_NilValue, + Rcpp::Nullable algo_parameters = R_NilValue) { typedef double NT; typedef Cartesian Kernel; @@ -50,26 +41,27 @@ Rcpp::List zono_approx (Rcpp::Reference Z, Rcpp::Nullable fit_ratio = R_Ni 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; + 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; NT ratio = std::numeric_limits::signaling_NaN(); - MT X(n, 2*k); + MT X(n, 2 * k); X << Rcpp::as(Z.field("G")).transpose(), -Rcpp::as(Z.field("G")).transpose(); - Eigen::JacobiSVD svd(X*X.transpose(), Eigen::ComputeFullU | Eigen::ComputeFullV); - MT G(k, 2*n); - G << Rcpp::as(Z.field("G"))*svd.matrixU(), Rcpp::as(Z.field("G"))*svd.matrixU(); + Eigen::JacobiSVD svd(X * X.transpose(), Eigen::ComputeFullU | Eigen::ComputeFullV); + MT G(k, 2 * n); + G << Rcpp::as(Z.field("G")) * svd.matrixU(), Rcpp::as(Z.field("G")) * svd.matrixU(); VT Gred_ii = G.transpose().cwiseAbs().rowwise().sum(); - MT A(n, 2*n); - A << -MT::Identity(n,n), MT::Identity(n,n); - MT Mat(2*n, n+1); + 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(); + Mat << Gred_ii, A.transpose() * svd.matrixU().transpose(); Hpolytope HP; - HP.init(n, A.transpose()*svd.matrixU().transpose(), Gred_ii); + HP.init(n, A.transpose() * svd.matrixU().transpose(), Gred_ii); if (fit_ratio.isNotNull() && Rcpp::as(fit_ratio)) { NT vol_red = std::abs(svd.matrixU().determinant()); @@ -77,58 +69,64 @@ Rcpp::List zono_approx (Rcpp::Reference Z, Rcpp::Nullable fit_ratio = R_Ni vol_red *= 2.0 * Gred_ii(i); } - if (error.isNotNull()) e = Rcpp::as(error); - if (walk_length.isNotNull()) walkL = Rcpp::as(walk_length); - if (rounding.isNotNull()) round = Rcpp::as(rounding); - if (!random_walk.isNotNull() || Rcpp::as(random_walk).compare(std::string("BiW")) == 0) { + walkL = (!Rcpp::as(algo_parameters).containsElementNamed("walk_length")) ? 1 : Rcpp::as( + Rcpp::as(algo_parameters)["walk_length"]); + e = (!Rcpp::as(algo_parameters).containsElementNamed("error")) ? 0.1 : Rcpp::as( + Rcpp::as(algo_parameters)["error"]); + round = (!Rcpp::as(algo_parameters).containsElementNamed("rounding")) ? false : Rcpp::as( + Rcpp::as(algo_parameters)["rounding"]); + //if (rounding.isNotNull()) round = Rcpp::as(rounding); + if (!Rcpp::as(algo_parameters).containsElementNamed("random_walk")) { billiard = true; win_len = 170; NN = 125; - } else if(Rcpp::as(random_walk).compare(std::string("RDHR")) == 0) { + } else if (Rcpp::as(Rcpp::as(algo_parameters)["random_walk"]).compare(std::string("BiW")) == 0) { + billiard = true; + win_len = 170; + NN = 125; + } else if (Rcpp::as(Rcpp::as(algo_parameters)["random_walk"]).compare(std::string("RDHR")) == 0) { rdhr = true; - } else if (Rcpp::as(random_walk).compare(std::string("CDHR")) == 0) { + } else if (Rcpp::as(Rcpp::as(algo_parameters)["random_walk"]).compare(std::string("CDHR")) == 0) { cdhr = true; - } else if (Rcpp::as(random_walk).compare(std::string("BaW")) == 0) { + } else if (Rcpp::as(Rcpp::as(algo_parameters)["random_walk"]).compare(std::string("Baw")) == 0) { ball_walk = true; - }else { + } else { throw Rcpp::exception("Unknown walk type!"); } - if(parameters.isNotNull()) { - if (Rcpp::as(parameters).containsElementNamed("BW_rad")) { - delta = Rcpp::as(Rcpp::as(parameters)["BW_rad"]); - } - if (Rcpp::as(parameters).containsElementNamed("Window")) { - win_len = Rcpp::as(Rcpp::as(parameters)["Window"]); - } - if (Rcpp::as(parameters).containsElementNamed("lb")) { - lb = Rcpp::as(Rcpp::as(parameters)["lb"]); - } - if (Rcpp::as(parameters).containsElementNamed("ub")) { - ub = Rcpp::as(Rcpp::as(parameters)["ub"]); - } - if (Rcpp::as(parameters).containsElementNamed("nu")) { - nu = Rcpp::as(Rcpp::as(parameters)["nu"]); - } - if (Rcpp::as(parameters).containsElementNamed("N")) { - NN = Rcpp::as(Rcpp::as(parameters)["N"]); - } - if (Rcpp::as(parameters).containsElementNamed("minmaxW")) { - win2 = Rcpp::as(Rcpp::as(parameters)["minmaxW"]); - } - if (Rcpp::as(parameters).containsElementNamed("prob")) { - p = Rcpp::as(Rcpp::as(parameters)["prob"]); - } - if (Rcpp::as(parameters).containsElementNamed("alpha")) { - alpha = Rcpp::as(Rcpp::as(parameters)["alpha"]); - } - if (Rcpp::as(parameters).containsElementNamed("hpoly")) { - hpoly = Rcpp::as(Rcpp::as(parameters)["hpoly"]); - } - if (Rcpp::as(parameters).containsElementNamed("diameter")) { - diam = Rcpp::as(Rcpp::as(parameters)["diameter"]); - } + if (Rcpp::as(algo_parameters).containsElementNamed("BW_rad")) { + delta = Rcpp::as(Rcpp::as(algo_parameters)["BW_rad"]); + } + if (Rcpp::as(algo_parameters).containsElementNamed("len_win")) { + win_len = Rcpp::as(Rcpp::as(algo_parameters)["len_win"]); + } + if (Rcpp::as(algo_parameters).containsElementNamed("lb")) { + lb = Rcpp::as(Rcpp::as(algo_parameters)["lb"]); + } + if (Rcpp::as(algo_parameters).containsElementNamed("ub")) { + ub = Rcpp::as(Rcpp::as(algo_parameters)["ub"]); + } + if (Rcpp::as(algo_parameters).containsElementNamed("nu")) { + nu = Rcpp::as(Rcpp::as(algo_parameters)["nu"]); + } + if (Rcpp::as(algo_parameters).containsElementNamed("N")) { + NN = Rcpp::as(Rcpp::as(algo_parameters)["N"]); + } + if (Rcpp::as(algo_parameters).containsElementNamed("minmaxW")) { + win2 = Rcpp::as(Rcpp::as(algo_parameters)["minmaxW"]); + } + if (Rcpp::as(algo_parameters).containsElementNamed("prob")) { + p = Rcpp::as(Rcpp::as(algo_parameters)["prob"]); + } + if (Rcpp::as(algo_parameters).containsElementNamed("alpha")) { + alpha = Rcpp::as(Rcpp::as(algo_parameters)["alpha"]); + } + if (Rcpp::as(algo_parameters).containsElementNamed("hpoly")) { + hpoly = Rcpp::as(Rcpp::as(algo_parameters)["hpoly"]); + } + if (Rcpp::as(algo_parameters).containsElementNamed("L")) { + diam = Rcpp::as(Rcpp::as(algo_parameters)["L"]); } zonotope ZP; @@ -136,35 +134,38 @@ Rcpp::List zono_approx (Rcpp::Reference Z, Rcpp::Nullable fit_ratio = R_Ni unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); // the random engine with this seed - typedef boost::mt19937 RNGType; + typedef boost::mt19937 RNGType; RNGType rng(seed); boost::random::uniform_real_distribution<>(urdist); - boost::random::uniform_real_distribution<> urdist1(-1,1); - - std::pair InnerB; - - if(inner_ball.isNotNull()) { //if it is given as an input - // store internal point hat is given as input - Rcpp::NumericVector InnerVec = Rcpp::as(inner_ball); - std::vector temp_p; - for (unsigned int j=0; j urdist1(-1, 1); + + + std::pair InnerB; + InnerB.second = -1.0; + if (Rcpp::as(algo_parameters).containsElementNamed("inner_ball")) { + if (Rcpp::as(Rcpp::as(algo_parameters)["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(algo_parameters)["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."); + } + } else { InnerB = ZP.ComputeInnerBall(); } - if (billiard && diam < 0.0){ + 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 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; @@ -172,12 +173,13 @@ Rcpp::List zono_approx (Rcpp::Reference Z, Rcpp::Nullable fit_ratio = R_Ni vol = vol_cooling_balls(ZP, var, var_ban, InnerB); } 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); + 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); } - ratio = std::pow(vol_red / vol, 1.0/NT(n)); + ratio = std::pow(vol_red / vol, 1.0 / NT(n)); } - return Rcpp::List::create(Rcpp::Named("Mat") = Rcpp::wrap(Mat) , Rcpp::Named("fit_ratio") = ratio); + return Rcpp::List::create(Rcpp::Named("Mat") = Rcpp::wrap(Mat), Rcpp::Named("fit_ratio") = ratio); } diff --git a/R-proj/tests/testthat/test_Hvol.R b/R-proj/tests/testthat/test_Hvol.R index 685ddd75d..8dafeb6ab 100644 --- a/R-proj/tests/testthat/test_Hvol.R +++ b/R-proj/tests/testthat/test_Hvol.R @@ -10,7 +10,7 @@ Hruntest <- function(P, name_string, exactvol, tol, num_of_exps, alg){ if (alg == "CB") { vol = vol + volume(P) } else { - vol = vol + volume(P, algo = "CG") + vol = vol + volume(P, algo = list("algorithm" = "CG")) } } vol = vol / num_of_exps diff --git a/R-proj/tests/testthat/test_Vvol.R b/R-proj/tests/testthat/test_Vvol.R index d29a04cf7..3427cf7ba 100644 --- a/R-proj/tests/testthat/test_Vvol.R +++ b/R-proj/tests/testthat/test_Vvol.R @@ -9,7 +9,7 @@ Vruntest <- function(P, name_string, exactvol, tol, num_of_exps, algorithm){ if (algorithm == "CB") { vol = vol + volume(P) } else { - vol = vol + volume(P, error=0.1, algo = "CG") + vol = vol + volume(P, algo = list("algorithm" = "CG", "error" = 0.1)) } } vol = vol / num_of_exps diff --git a/R-proj/tests/testthat/test_Zvol.R b/R-proj/tests/testthat/test_Zvol.R index 18306e71e..8da513ec1 100644 --- a/R-proj/tests/testthat/test_Zvol.R +++ b/R-proj/tests/testthat/test_Zvol.R @@ -8,9 +8,9 @@ Zruntest <- function(P, name_string, tol, num_of_exps, algo){ vol = 0 for (j in 1:num_of_exps) { if (algo == "CB") { - vol = vol + volume(P, rounding=TRUE) + vol = vol + volume(P, rounding=FALSE) } else { - vol = vol + volume(P, error=0.1, algo = "CG", rounding=TRUE) + vol = vol + volume(P, algo = list("algorithm" = "CG", "error" = 0.1), rounding=TRUE) } } vol = vol / num_of_exps diff --git a/R-proj/tests/testthat/test_sampling.R b/R-proj/tests/testthat/test_sampling.R index e626cf2ff..6df5df4d9 100644 --- a/R-proj/tests/testthat/test_sampling.R +++ b/R-proj/tests/testthat/test_sampling.R @@ -4,9 +4,9 @@ library(volesti) runsample <- function(P, name_string, dist){ if (dist == "uniform") { - p = sample_points(P) + p = sample_points(P, n = 100) } else { - p = sample_points(P, distribution = "gaussian") + p = sample_points(P, n = 100, distribution = list("density" = "gaussian")) } if (length(p[is.nan(p)])>0 | length(p[is.infinite(p)])>0) { res = 0 diff --git a/cran_gen/NEWS.md b/cran_gen/NEWS.md index 2a5eaafc1..872400077 100644 --- a/cran_gen/NEWS.md +++ b/cran_gen/NEWS.md @@ -19,5 +19,7 @@ * New volume computation algorithm. * Billiard walk for uniform sampling. * Modified exact volume computation function. +* Implementation and evaluation of PCA method for zonotope approximation. +* Boundary sampling. * Improved functionality for finance applications. * Improved names for functions and input variables. diff --git a/include/generators/h_polytopes_gen.h b/include/generators/h_polytopes_gen.h index 047d77c83..f0e273d64 100644 --- a/include/generators/h_polytopes_gen.h +++ b/include/generators/h_polytopes_gen.h @@ -12,15 +12,19 @@ #include "samplers.h" template -Polytope random_hpoly(unsigned int dim, unsigned int m) { +Polytope random_hpoly(unsigned int dim, unsigned int m, double seed = std::numeric_limits::signaling_NaN()) { typedef typename Polytope::MT MT; typedef typename Polytope::VT VT; typedef typename Polytope::NT NT; typedef typename Polytope::PolytopePoint Point; - 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)) { + unsigned rng_seed = seed; + rng.seed(rng_seed); + } boost::random::uniform_real_distribution<> urdist1(-10, 10); Point p(dim); typename std::vector::iterator pit; diff --git a/include/generators/v_polytopes_gen.h b/include/generators/v_polytopes_gen.h index 73e148bbf..10191b647 100644 --- a/include/generators/v_polytopes_gen.h +++ b/include/generators/v_polytopes_gen.h @@ -25,27 +25,44 @@ void removeRow(MT &matrix, unsigned int rowToRemove) } template -Polytope random_vpoly(unsigned int dim, unsigned int k) { +Polytope random_vpoly(unsigned int dim, unsigned int k, double seed = std::numeric_limits::signaling_NaN()) { typedef typename Polytope::MT MT; typedef typename Polytope::VT VT; typedef typename Polytope::NT NT; typedef typename Polytope::PolytopePoint Point; - 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)) { + unsigned rng_seed = seed; + rng.seed(rng_seed); + } + boost::normal_distribution<> rdist(0,1); - Point p; typename std::vector::iterator pit; MT V(k, dim); unsigned int j; + + std::vector Xs(dim,0); + NT normal = NT(0); + for (unsigned int i = 0; i < k; ++i) { - p = get_direction(dim); - pit = p.iter_begin(); - j = 0; - for ( ; pit!=p.iter_end(); ++pit, ++j) { - V(i,j) = *pit; + + normal = NT(0); + for (unsigned int i=0; i -Polytope random_vpoly_incube(unsigned int d, unsigned int k) { +Polytope random_vpoly_incube(unsigned int d, unsigned int k, double seed = std::numeric_limits::signaling_NaN()) { typedef typename Polytope::MT MT; typedef typename Polytope::VT VT; @@ -72,8 +89,12 @@ Polytope random_vpoly_incube(unsigned int d, unsigned int k) { conv_mem = (REAL *) malloc(k * sizeof(*conv_mem)); colno_mem = (int *) malloc(k * sizeof(*colno_mem)); - 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)) { + unsigned rng_seed = seed; + rng.seed(rng_seed); + } boost::random::uniform_real_distribution<> urdist1(-1, 1); Point p(d); diff --git a/include/generators/z_polytopes_gen.h b/include/generators/z_polytopes_gen.h index 3aa818df0..4c9b3dedd 100644 --- a/include/generators/z_polytopes_gen.h +++ b/include/generators/z_polytopes_gen.h @@ -12,14 +12,18 @@ #include "samplers.h" template -Polytope gen_zonotope_gaussian(int dim, int m) { +Polytope gen_zonotope_gaussian(int dim, int m, double seed = std::numeric_limits::signaling_NaN()) { typedef typename Polytope::MT MT; typedef typename Polytope::VT VT; typedef typename Polytope::NT NT; - 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)) { + unsigned rng_seed = seed; + rng.seed(rng_seed); + } boost::normal_distribution<> rdist(0, 1); boost::normal_distribution<> rdist2(50, 33.3); @@ -51,14 +55,18 @@ Polytope gen_zonotope_gaussian(int dim, int m) { template -Polytope gen_zonotope_uniform(int dim, int m) { +Polytope gen_zonotope_uniform(int dim, int m, double seed = std::numeric_limits::signaling_NaN()) { typedef typename Polytope::MT MT; typedef typename Polytope::VT VT; typedef typename Polytope::NT NT; - 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)) { + unsigned rng_seed = seed; + rng.seed(rng_seed); + } boost::normal_distribution<> rdist(0, 1); boost::random::uniform_real_distribution<> urdist1(0, 100); @@ -83,14 +91,18 @@ Polytope gen_zonotope_uniform(int dim, int m) { template -Polytope gen_zonotope_exponential(int dim, int m) { +Polytope gen_zonotope_exponential(int dim, int m, double seed = std::numeric_limits::signaling_NaN()) { typedef typename Polytope::MT MT; typedef typename Polytope::VT VT; typedef typename Polytope::NT NT; - 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)) { + unsigned rng_seed = seed; + rng.seed(rng_seed); + } boost::normal_distribution<> rdist(0, 1); boost::normal_distribution<> expdist(1.0/30.0); diff --git a/include/rounding.h b/include/rounding.h index 65d6a6389..90e5b8f5d 100644 --- a/include/rounding.h +++ b/include/rounding.h @@ -5,6 +5,10 @@ //Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 program. +// The functions in this header file call Bojan Nikolic implementation +// of Todd and Yildirim algorithm in "On Khachiyan's Algorithm for the Computation of Minimum Volume Enclosing Ellipsoids", 2005 +// Initial version 2010, part of BNMin1 and is licensed under GNU General Public License version 2 + // Licensed under GNU LGPL.3, see LICENCE file #ifndef ROUNDING_H @@ -30,27 +34,19 @@ std::pair rounding_min_ellipsoid(Polytope &P , const std::pair20*domension // 2. Generate the first random point in P // Perform random walk on random point in the Chebychev ball - Point p = get_point_on_Dsphere(n, radius); + Point p = get_point_in_Dsphere(n, radius); p = p + c; //use a large walk length e.g. 1000 - Parameters var2 = var; - var2.cdhr_walk = true; - var2.ball_walk = var2.rdhr_walk = var2.bill_walk = false; - rand_point_generator(P, p, 1, 50*n, randPoints, var2); + rand_point_generator(P, p, 1, 10*n, randPoints, var); // 3. Sample points from P unsigned int num_of_samples = 10*n;//this is the number of sample points will used to compute min_ellipoid randPoints.clear(); - rand_point_generator(P, p, num_of_samples, 10 + n/10, randPoints, var2); - /*NT current_dist, max_dist; - for(typename std::list::iterator pit=randPoints.begin(); pit!=randPoints.end(); ++pit){ - current_dist=(*pit-c).squared_length(); - if(current_dist>max_dist){ - max_dist=current_dist; - } + if (var.bill_walk) { + rand_point_generator(P, p, num_of_samples, 5, randPoints, var); + } else { + rand_point_generator(P, p, num_of_samples, 10 + n / 10, randPoints, var); } - max_dist=std::sqrt(max_dist); - R=max_dist/radius;*/ } // Store points in a matrix to call Khachiyan algorithm for the minimum volume enclosing ellipsoid diff --git a/include/samplers/sample_only.h b/include/samplers/sample_only.h index 57905f315..7bbcacd00 100644 --- a/include/samplers/sample_only.h +++ b/include/samplers/sample_only.h @@ -24,18 +24,16 @@ template void sampling_only(PointList &randPoints, Polytope &P, const unsigned int walk_len, - const unsigned int rnum, bool gaussian, const NT &a, const Point &internal_point, - UParameters const& var1, GParameters const& var2) { + const unsigned int rnum, bool gaussian, const NT &a, const bool boundary, + const Point &internal_point, UParameters const& var1, GParameters const& var2) { typedef typename UParameters::RNGType RNGType; - unsigned int n = var1.n; Point p = internal_point; - Point q = get_point_on_Dsphere(n, var1.che_rad); - p=p+q; - rand_point_generator(P, p, 1, 50 * n, randPoints, var1); randPoints.clear(); - if (!gaussian){ + 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); diff --git a/test/cooling_bodies_bill_test.cpp b/test/cooling_bodies_bill_test.cpp index 895dec6e1..6708b680d 100644 --- a/test/cooling_bodies_bill_test.cpp +++ b/test/cooling_bodies_bill_test.cpp @@ -59,11 +59,16 @@ void test_cool_bodies(Polytope &HP, NT expected, NT tolerance=0.1, bool round = std::pair InnerBall; if (round) { + HP.normalize(); InnerBall = HP.ComputeInnerBall(); - vars var2(rnum,n,walk_len,n_threads,err,e,0,0,0,0,0.0,rng, + HP.comp_diam(diameter, InnerBall.second); + diameter*=2.0; + + vars var2(rnum,n,walk_len,n_threads,err,e,0,0,0,0,diameter,rng, urdist,urdist1,-1.0,false,false,false,false,false,false,false,false,true); std::pair res_round = rounding_min_ellipsoid(HP , InnerBall, var2); round_val = res_round.first; + diameter = -1.0; } InnerBall = HP.ComputeInnerBall(); diff --git a/test/vol.cpp b/test/vol.cpp index 22cb62ac3..c8efc5e66 100644 --- a/test/vol.cpp +++ b/test/vol.cpp @@ -89,7 +89,8 @@ int main(const int argc, const char** argv) gaussian_sam = false, hpoly = false, billiard=false, - win2 = false; + win2 = false, + boundary = false; //this is our polytope Hpolytope HP; @@ -633,11 +634,11 @@ int main(const int argc, const char** argv) double tstart11 = (double)clock()/(double)CLOCKS_PER_SEC; if (Zono) { - sampling_only(randPoints, ZP, walk_len, nsam, gaussian_sam, a, InnerBall.first, var1, var2); + sampling_only(randPoints, ZP, walk_len, nsam, gaussian_sam, a, boundary, InnerBall.first, var1, var2); } else if (!Vpoly) { - sampling_only(randPoints, HP, walk_len, nsam, gaussian_sam, a, InnerBall.first, var1, var2); + sampling_only(randPoints, HP, walk_len, nsam, gaussian_sam, a, boundary, InnerBall.first, var1, var2); } else { - sampling_only(randPoints, VP, walk_len, nsam, gaussian_sam, a, InnerBall.first, var1, var2); + sampling_only(randPoints, VP, walk_len, nsam, gaussian_sam, a, boundary, InnerBall.first, var1, var2); } double tstop11 = (double)clock()/(double)CLOCKS_PER_SEC; if(verbose) std::cout << "Sampling time: " << tstop11 - tstart11 << std::endl;