diff --git a/R-proj/DESCRIPTION b/R-proj/DESCRIPTION index ffda1a9fa..3ee4bdfc3 100644 --- a/R-proj/DESCRIPTION +++ b/R-proj/DESCRIPTION @@ -12,7 +12,7 @@ Description: Provides an R interface for 'volesti' C++ package. 'volesti' comput for sampling, rounding and rotating polytopes. Moreover, 'volesti' provides algorithms for estimating copulas useful in computational finance. Version: 1.1.0 -Date: 2020-01-23 +Date: 2020-05-29 Maintainer: Vissarion Fisikopoulos Depends: Rcpp (>= 0.12.17) Imports: methods, stats diff --git a/R-proj/NAMESPACE b/R-proj/NAMESPACE index ee585568c..acfb34451 100644 --- a/R-proj/NAMESPACE +++ b/R-proj/NAMESPACE @@ -22,7 +22,7 @@ export(volume) export(zonotope_approximation) exportPattern("^[[:alpha:]]+") importFrom("methods","new") -importFrom("stats", "cov") +importFrom("stats","cov") importFrom("utils","read.csv") importFrom(Rcpp,evalCpp) importFrom(Rcpp,loadModule) diff --git a/R-proj/R/RcppExports.R b/R-proj/R/RcppExports.R index 1362b3ecc..6250dce56 100644 --- a/R-proj/R/RcppExports.R +++ b/R-proj/R/RcppExports.R @@ -3,10 +3,11 @@ #' Construct a copula using uniform sampling from the unit simplex #' -#' Given two families of parallel hyperplanes or a family of parallel hyperplanes and a family of concentric ellispoids centered at the origin intersecting the canonical simplex, this function uniformly samples from the canonical simplex and construct an approximation of the bivariate probability distribution, called copula. +#' Given two families of parallel hyperplanes or a family of parallel hyperplanes and a family of concentric ellispoids centered at the origin intersecting the canonical simplex, this function uniformly samples from the canonical simplex and construct an approximation of the bivariate probability distribution, called copula (see \url{https://en.wikipedia.org/wiki/Copula_(probability_theory)}). +#' At least two families of hyperplanes or one family of hyperplanes and one family of ellipsoids have to be given as input. #' -#' @param r1 A \eqn{d}-dimensional vector that describes the direction of the first family of parallel hyperplanes. -#' @param r2 Optional. A \eqn{d}-dimensional vector that describes the direction of the second family of parallel hyperplanes. +#' @param r1 The \eqn{d}-dimensional normal vector of the first family of parallel hyperplanes. +#' @param r2 Optional. The \eqn{d}-dimensional normal vector of the second family of parallel hyperplanes. #' @param sigma Optional. The \eqn{d\times d} symmetric positive semidefine matrix that describes the family of concentric ellipsoids centered at the origin. #' @param m The number of the slices for the copula. The default value is 100. #' @param n The number of points to sample. The default value is \eqn{5\cdot 10^5}. @@ -15,7 +16,7 @@ #' @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.} #' -#' @return A \eqn{numSlices\times numSlices} numerical matrix that corresponds to a copula. +#' @return A \eqn{m\times m} numerical matrix that corresponds to a copula. #' @examples #' # compute a copula for two random families of parallel hyperplanes #' h1 = runif(n = 10, min = 1, max = 1000) @@ -32,7 +33,7 @@ #' cop = copula(r1 = h, sigma = E, m = 10, n = 100000) #' #' @export -copula <- function(r1 = NULL, r2 = NULL, sigma = NULL, m = NULL, n = NULL, seed = NULL) { +copula <- function(r1, r2 = NULL, sigma = NULL, m = NULL, n = NULL, seed = NULL) { .Call(`_volesti_copula`, r1, r2, sigma, m, n, seed) } @@ -59,7 +60,7 @@ copula <- function(r1 = NULL, r2 = NULL, sigma = NULL, m = NULL, n = NULL, seed #' # 100 uniform points from the 2-d unit ball #' points = direct_sampling(n = 100, body = list("type" = "ball", "dimension" = 2)) #' @export -direct_sampling <- function(body = NULL, n = NULL, seed = NULL) { +direct_sampling <- function(body, n, seed = NULL) { .Call(`_volesti_direct_sampling`, body, n, seed) } @@ -70,11 +71,14 @@ direct_sampling <- function(body = NULL, n = NULL, seed = NULL) { #' #' @param P A polytope #' +#' @references \cite{E. Gover and N. Krikorian, +#' \dQuote{Determinants and the Volumes of Parallelotopes and Zonotopes,} \emph{Linear Algebra and its Applications, 433(1), 28 - 40,} 2010.} +#' #' @return The exact volume of the input polytope, for zonotopes, simplices in V-representation and polytopes with known exact volume #' @examples #' #' # compute the exact volume of a 5-dimensional zonotope defined by the Minkowski sum of 10 segments -#' Z = gen_rand_zonotope(5, 10) +#' Z = gen_rand_zonotope(2, 5) #' vol = exact_vol(Z) #' #' \donttest{# compute the exact volume of a 2-d arbitrary simplex @@ -91,9 +95,9 @@ exact_vol <- function(P) { .Call(`_volesti_exact_vol`, P) } -#' Compute the percentage of the volume of the unit simplex that is contained in the intersection of a half-space and the unit simplex. +#' Compute the percentage of the volume of the simplex that is contained in the intersection of a half-space and the simplex. #' -#' A half-space \eqn{H} is given as a pair of a vector \eqn{a\in R^d} and a scalar \eqn{z0\in R} s.t.: \eqn{a^Tx\leq z0}. This function calls the Ali's version of the Varsi formula to compute a frustum of the unit simplex. +#' A half-space \eqn{H} is given as a pair of a vector \eqn{a\in R^d} and a scalar \eqn{z0\in R} s.t.: \eqn{a^Tx\leq z0}. This function calls the Ali's version of the Varsi formula to compute a frustum of the simplex. #' #' @param a A \eqn{d}-dimensional vector that defines the direction of the hyperplane. #' @param z0 The scalar that defines the half-space. @@ -104,7 +108,7 @@ exact_vol <- function(P) { #' @references \cite{Ali, Mir M., #' \dQuote{Content of the frustum of a simplex,} \emph{ Pacific J. Math. 48, no. 2, 313--322,} 1973.} #' -#' @return The percentage of the volume of the unit simplex that is contained in the intersection of a given half-space and the unit simplex. +#' @return The percentage of the volume of the simplex that is contained in the intersection of a given half-space and the simplex. #' #' @examples #' # compute the frustum of H: -x1+x2<=0 @@ -118,13 +122,12 @@ frustum_of_simplex <- function(a, z0) { #' Compute an inscribed ball of a convex polytope #' -#' For a H-polytope described by a \eqn{m\times d} matrix \eqn{A} and a \eqn{m}-dimensional vector \eqn{b}, s.t.: \eqn{Ax\leq b}, this function computes the largest inscribed ball (Chebychev ball) by solving the corresponding linear program. -#' For a V-polytope \eqn{d+1} vertices, that define a full dimensional simplex, picked at random and the largest inscribed ball of the simplex is computed. -#' For a zonotope \eqn{P} we compute the minimum \eqn{r} s.t.: \eqn{ r e_i \in P} for all \eqn{i=1, \dots ,d}. Then the ball centered at the origin with radius \eqn{r/ \sqrt{d}} is an inscribed ball. +#' For a H-polytope described by a \eqn{m\times d} matrix \eqn{A} and a \eqn{m}-dimensional vector \eqn{b}, s.t.: \eqn{P=\{x\ |\ Ax\leq b\} }, this function computes the largest inscribed ball (Chebychev ball) by solving the corresponding linear program. +#' For both zonotopes and V-polytopes the function computes the minimum \eqn{r} s.t.: \eqn{ r e_i \in P} 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 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 or (b) Vpolytope or (c) Zonotope or (d) VpolytopeIntersection. #' -#' @return A \eqn{d+1}-dimensional vector that describes the inscribed ball. The first \eqn{d} coordinates corresponds to the center of the ball and the last one to the radius. +#' @return A \eqn{(d+1)}-dimensional vector that describes the inscribed ball. The first \eqn{d} coordinates corresponds to the center of the ball and the last one to the radius. #' #' @examples #' # compute the Chebychev ball of the 2d unit simplex @@ -148,8 +151,7 @@ inner_ball <- function(P) { #' @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 Optional. A fixed seed for the random polytope generator. #' -#' @section warning: -#' Do not use this function. +#' @keywords internal #' #' @return A numerical matrix describing the requested polytope poly_gen <- function(kind_gen, Vpoly_gen, Zono_gen, dim_gen, m_gen, seed = NULL) { @@ -162,8 +164,7 @@ poly_gen <- function(kind_gen, Vpoly_gen, Zono_gen, dim_gen, m_gen, seed = NULL) #' @param T Optional. A rotation matrix. #' @param seed Optional. A fixed seed for the random linear map generator. #' -#' @section warning: -#' Do not use this function. +#' @keywords internal #' #' @return A matrix that describes the rotated polytope rotating <- function(P, T = NULL, seed = NULL) { @@ -175,41 +176,40 @@ rotating <- function(P, T = NULL, seed = NULL) { #' @param P A convex polytope (H- or V-representation or zonotope). #' @param seed Optional. A fixed seed for the number generator. #' -#' @section warning: -#' Do not use this function. +#' @keywords internal #' -#' @return A numerical matrix that describes the rounded polytope and contains the round value. +#' @return A numerical matrix that describes the rounded polytope, a numerical matrix of the inverse linear transofmation that is applied on the input polytope, the numerical vector the the input polytope is shifted and the determinant of the matrix of the linear transformation that is applied on the input polytope. rounding <- function(P, seed = NULL) { .Call(`_volesti_rounding`, P, seed) } -#' Sample uniformly or normally distributed points from a convex Polytope (H-polytope, V-polytope or a zonotope). +#' Sample uniformly or normally distributed points from a convex Polytope (H-polytope, V-polytope, zonotope or intersection of two V-polytopes). #' -#' Sample n points with uniform or multidimensional spherical gaussian -with a mode at any point- target distribution. +#' Sample n points with uniform or multidimensional spherical gaussian -with a mode at any point- as the target distribution. #' -#' @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 P A convex polytope. It is an object from class (a) Hpolytope or (b) Vpolytope or (c) Zonotope or (d) VpolytopeIntersection. #' @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{walk_length} }{ The number of the steps per generated point for the random walk. The default value is \eqn{5} for \code{'BiW'} and \eqn{\lfloor 10 + d/10\rfloor} otherwise.} #' \item{\code{nburns} }{ The number of points to burn before start sampling.} -#' \item{\code{starting_point} }{ A \eqn{d}-dimensional numerical vector that declares a starting point in the interior of the polytope for the random walk. The default choice is the center of the Chebychev ball.} +#' \item{\code{starting_point} }{ A \eqn{d}-dimensional numerical vector that declares a starting point in the interior of the polytope for the random walk. The default choice is the center of the ball as that one computed by the function \code{inner_ball()}.} #' \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{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{mode} }{ A \eqn{d}-dimensional numerical vector that declares the mode of the Gaussian distribution. The default choice is the center of the Chebychev ball.} +#' \item{\code{mode} }{ A \eqn{d}-dimensional numerical vector that declares the mode of the Gaussian distribution. The default choice is the center of the as that one computed by the function \code{inner_ball()}.} #' } #' @param seed Optional. A fixed seed for the number generator. #' #' @return A \eqn{d\times n} matrix that contains, column-wise, the sampled points from the convex polytope P. #' @examples -#' # uniform distribution from the 3d unit cube in V-representation using ball walk -#' P = gen_cube(3, 'V') +#' # uniform distribution from the 3d unit cube in H-representation using ball walk +#' P = gen_cube(3, 'H') #' 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 @@ -220,32 +220,32 @@ rounding <- function(P, seed = NULL) { #' #' # 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")) +#' points = sample_points(P, n = 100, random_walk = list("walk" = "BRDHR")) #' #' @export -sample_points <- function(P = NULL, n = NULL, random_walk = NULL, distribution = NULL, seed = NULL) { +sample_points <- function(P, n, random_walk = NULL, distribution = NULL, seed = NULL) { .Call(`_volesti_sample_points`, P, n, random_walk, distribution, seed) } -#' The main function for volume approximation of a convex Polytope (H-polytope, V-polytope or a zonotope) +#' The main function for volume approximation of a convex Polytope (H-polytope, V-polytope, zonotope or intersection of two V-polytopes) #' -#' 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. +#' For the volume approximation can be used three algorithms. Either CoolingBodies (CB) or SequenceOfBalls (SOB) or CoolingGaussian (CG). An 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{P=\{x\ |\ 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 P A convex polytope. It is an object from class a) Hpolytope or b) Vpolytope or c) Zonotope or d) VpolytopeIntersection. #' @param settings Optional. A list that declares which algorithm, random walk and values of parameters to use, as follows: #' \itemize{ -#' \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{algorithm} }{ A string to set the algorithm to use: a) \code{'CB'} for CB algorithm, b) \code{'SoB'} for SOB algorithm or b) \code{'CG'} for CG algorithm. 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 SOB algorithm 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. For CB and SOB algorithms the default walk is \code{'CDHR'} for H-polytopes and \code{'BiW'} for the other representations. For CG algorithm the default walk is \code{'CDHR'} for H-polytopes and \code{'RDHR'} 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{win_len} }{ The length of the sliding window for CG algorithm. The default value is \eqn{500+4dimension^2}.} -#' \item{\code{hpoly} }{ A boolean parameter to use H-polytopes in MMC of CB algorithm. The default value is \code{FALSE}.} +#' \item{\code{win_len} }{ The length of the sliding window for CB or CG algorithm. The default value is \eqn{400+3d^2} for CB or \eqn{500+4d^2} for CG.} +#' \item{\code{hpoly} }{ A boolean parameter to use H-polytopes in MMC of CB algorithm when the input polytope is a zonotope. The default value is \code{TRUE} when the order of the zonotope is \eqn{<5}, otherwise it is \code{FALSE}.} #' } #' @param rounding Optional. A boolean parameter for rounding. The default value is \code{TRUE} for V-polytopes and \code{FALSE} otherwise. #' @param seed Optional. A fixed seed for the number generator. #' #' @references \cite{I.Z.Emiris and V. Fisikopoulos, -#' \dQuote{Practical polytope volume approximation,} \emph{ACM Trans. Math. Soft.,} 2014.}, +#' \dQuote{Practical polytope volume approximation,} \emph{ACM Trans. Math. Soft.,} 2018.}, #' @references \cite{A. Chalkis and I.Z.Emiris and V. Fisikopoulos, #' \dQuote{Practical Volume Estimation by a New Annealing Schedule for Cooling Convex Bodies,} \emph{CoRR, abs/1905.05494,} 2019.}, #' @references \cite{B. Cousins and S. Vempala, \dQuote{A practical volume algorithm,} \emph{Springer-Verlag Berlin Heidelberg and The Mathematical Programming Society,} 2015.} @@ -253,17 +253,19 @@ sample_points <- function(P = NULL, n = NULL, random_walk = NULL, distribution = #' #' @return The approximation of the volume of a convex polytope. #' @examples -#' # calling SOB algorithm for a H-polytope (2d unit simplex) -#' P = gen_simplex(2,'H') -#' vol = volume(P) #' -#' # calling CG algorithm for a V-polytope (3d simplex) -#' P = gen_simplex(2,'V') -#' vol = volume(P, settings = list("algorithm" = "CG")) +#' # calling SOB algorithm for a H-polytope (3d unit simplex) +#' HP = gen_cube(3,'H') +#' vol = volume(HP) +#' +#' # calling CG algorithm for a V-polytope (2d simplex) +#' VP = gen_simplex(2,'V') +#' vol = volume(VP, settings = 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, settings = list("random_walk" = "RDHR", "walk_length" = 5)) +#' vol = volume(Z, settings = list("random_walk" = "RDHR", "walk_length" = 2)) +#' #' @export volume <- function(P, settings = NULL, rounding = NULL, seed = NULL) { .Call(`_volesti_volume`, P, settings, rounding, seed) @@ -276,8 +278,7 @@ volume <- function(P, settings = NULL, rounding = NULL, seed = NULL) { #' @param settings Optional. A list that declares the values of the parameters of CB algorithm. #' @param seed Optional. A fixed seed for the number generator. #' -#' @section warning: -#' Do not use this function. +#' @keywords internal #' #' @return A List that contains a numerical matrix that describes the PCA approximation as a H-polytope and the ratio of fitness. zono_approx <- function(Z, fit_ratio = NULL, settings = NULL, seed = NULL) { diff --git a/R-proj/R/compute_indicators.R b/R-proj/R/compute_indicators.R index 00d223693..fe82cdee3 100644 --- a/R-proj/R/compute_indicators.R +++ b/R-proj/R/compute_indicators.R @@ -1,8 +1,8 @@ #' Compute an indicator for each time period that describes the state of a market. #' -#' Given a matrix that contains row-wise the assets' returns and a sliding window W, this function computes an approximation of the joint distribution (copula) between portfolios' return and volatility in each time period implied by W. -#' For each copula it computes an indicator: large value corresponds to a crisis period and a small value to a normal period. -#' The periods over which the indicator is greater than 1 for more than 60 consecutives sliding windows are warnings and for more than 100 are crisis. The sliding window is shifted by one day. +#' Given a matrix that contains row-wise the assets' returns and a sliding window \code{win_length}, this function computes an approximation of the joint distribution (copula, e.g. see \url{https://en.wikipedia.org/wiki/Copula_(probability_theory)}) between portfolios' return and volatility in each time period defined by \code{win_len}. +#' For each copula it computes an indicator: If the indicator is large it corresponds to a crisis period and if it is small it corresponds to a normal period. +#' In particular, the periods over which the indicator is greater than 1 for more than 60 consecutive sliding windows are warnings and for more than 100 are crisis. The sliding window is shifted by one day. #' #' @param returns A \eqn{d}-dimensional vector that describes the direction of the first family of parallel hyperplanes. #' @param win_length Optional. The length of the sliding window. The default value is 60. @@ -17,6 +17,11 @@ #' #' @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. #' +#' @examples +#' # simple example on random asset returns +#' asset_returns = replicate(10, rnorm(14)) +#' market_states_and_indicators = compute_indicators(asset_returns, 10, 10, 10000, 2, 3) +#' #' @export compute_indicators <- function(returns, win_length = NULL, m = NULL, n = NULL, nwarning = NULL, ncrisis = NULL, seed = NULL) { @@ -77,15 +82,22 @@ compute_indicators <- function(returns, win_length = NULL, m = NULL, n = NULL, n set_index = TRUE } else if (indicators[i]<1) { if(set_index){ - if(i-index+1 > nwarning && i-index+1 < ncrisis){ - col[index:i] = "warning" - } else if(i-index+1 > ncrisis) { - col[index:i] = "crisis" + if(i-index > nwarning-1 && i-index <= ncrisis-1){ + col[index:(i-1)] = "warning" + } else if(i-index > ncrisis-1) { + col[index:(i-1)] = "crisis" } } set_index = FALSE } } + if(set_index){ + if(N-index+1 > nwarning-1 && N-index+1 <= ncrisis-1){ + col[index:i] = "warning" + } else if(N-index+1 > ncrisis-1) { + col[index:i] = "crisis" + } + } return(list("indicators" = indicators, market_states = col)) diff --git a/R-proj/R/file_to_polytope.R b/R-proj/R/file_to_polytope.R index f2decb056..52f469d87 100644 --- a/R-proj/R/file_to_polytope.R +++ b/R-proj/R/file_to_polytope.R @@ -1,6 +1,7 @@ #' function to get an ine or an ext file and returns the corresponding polytope #' -#' For an ine file it generates the corresponding H-polytope. For an ext file it generates the corresponding V-polytope or zonotope. +#' For an ".ine" file it generates the corresponding H-polytope. For an ".ext" file it generates the corresponding V-polytope or zonotope. +#' For more details on those file formats see \url{https://github.com/GeomScale/volume_approximation/blob/develop/doc/cpp_interface.md#polytope-input}. #' #' @param path A string that containes the path to an ine or a ext file. The ine file desrcibes a H-polytope and ext file describes a V-polytope or a zonotope. #' @param zonotope A boolean parameter. It has to be TRUE when the path leads to an .ext file that describes a zonotope. @@ -16,6 +17,7 @@ #' @importFrom Rcpp evalCpp #' @importFrom Rcpp loadModule #' @importFrom "utils" "read.csv" +#' @importFrom "stats" "cov" #' @importFrom "methods" "new" #' @exportPattern "^[[:alpha:]]+" file_to_polytope <- function(path, zonotope){ diff --git a/R-proj/R/gen_cross.R b/R-proj/R/gen_cross.R index 5eceac232..0ab388f62 100644 --- a/R-proj/R/gen_cross.R +++ b/R-proj/R/gen_cross.R @@ -1,6 +1,6 @@ #' Generator function for cross polytopes #' -#' This function can be used to generate the \eqn{d}-dimensional cross polytope in H- or V-representation. +#' This function generates the \eqn{d}-dimensional cross polytope in H- or V-representation. #' #' @param dimension The dimension of the cross polytope. #' @param representation A string to declare the representation. It has to be \code{'H'} for H-representation or \code{'V'} for V-representation. @@ -8,7 +8,7 @@ #' @return A polytope class representing a cross polytope in H- or V-representation. #' @examples #' # generate a 10-dimensional cross polytope in H-representation -#' P = gen_cross(10, 'H') +#' P = gen_cross(5, 'H') #' #' # generate a 15-dimension cross polytope in V-representation #' P = gen_cross(15, 'V') diff --git a/R-proj/R/gen_cube.R b/R-proj/R/gen_cube.R index 581131256..bcbf1f025 100644 --- a/R-proj/R/gen_cube.R +++ b/R-proj/R/gen_cube.R @@ -1,6 +1,6 @@ #' Generator function for hypercubes #' -#' This function can be used to generate the \eqn{d}-dimensional unit hypercube \eqn{[-1,1]^d} in H- or V-representation. +#' This function generates the \eqn{d}-dimensional unit hypercube \eqn{[-1,1]^d} in H- or V-representation. #' #' @param dimension The dimension of the hypercube #' @param representation A string to declare the representation. It has to be \code{'H'} for H-representation or \code{'V'} for V-representation. @@ -11,7 +11,7 @@ #' P = gen_cube(10, 'H') #' #' # generate a 15-dimension hypercube in V-representation -#' P = gen_cube(15, 'V') +#' P = gen_cube(5, 'V') #' @export gen_cube <- function(dimension, representation) { diff --git a/R-proj/R/gen_prod_simplex.R b/R-proj/R/gen_prod_simplex.R index 2eda50bd9..a0505880c 100644 --- a/R-proj/R/gen_prod_simplex.R +++ b/R-proj/R/gen_prod_simplex.R @@ -1,6 +1,6 @@ #' Generator function for product of simplices #' -#' This function can be used to generate a \eqn{2d}-dimensional polytope that is defined as the product of two \eqn{d}-dimensional unit simplices in H-representation. +#' This function generates a \eqn{2d}-dimensional polytope that is defined as the product of two \eqn{d}-dimensional unit simplices in H-representation. #' #' @param dimension The dimension of the simplices. #' diff --git a/R-proj/R/gen_rand_hpoly.R b/R-proj/R/gen_rand_hpoly.R index 354e9d268..2785a28b4 100644 --- a/R-proj/R/gen_rand_hpoly.R +++ b/R-proj/R/gen_rand_hpoly.R @@ -1,6 +1,6 @@ #' Generator function for random H-polytopes #' -#' This function can be used to generate a \eqn{d}-dimensional polytope in H-representation with \eqn{m} facets. We pick \eqn{m} random hyperplanes tangent on the \eqn{d}-dimensional unit hypersphere as facets. +#' This function generates a \eqn{d}-dimensional polytope in H-representation with \eqn{m} facets. We pick \eqn{m} random hyperplanes tangent on the \eqn{d}-dimensional unit hypersphere as facets. #' #' @param dimension The dimension of the convex polytope. #' @param nfacets The number of the facets. diff --git a/R-proj/R/gen_rand_vpoly.R b/R-proj/R/gen_rand_vpoly.R index 3711938be..31cc561c0 100644 --- a/R-proj/R/gen_rand_vpoly.R +++ b/R-proj/R/gen_rand_vpoly.R @@ -1,6 +1,6 @@ #' Generator function for random V-polytopes #' -#' This function can be used to generate a \eqn{d}-dimensional polytope in V-representation with \eqn{m} vertices. We pick \eqn{m} random points from the boundary of the \eqn{d}-dimensional unit hypersphere as vertices. +#' This function generates a \eqn{d}-dimensional polytope in V-representation with \eqn{m} vertices. We pick \eqn{m} random points from the boundary of the \eqn{d}-dimensional unit hypersphere as vertices. #' #' @param dimension The dimension of the convex polytope. #' @param nvertices The number of the vertices. diff --git a/R-proj/R/gen_rand_zonotope.R b/R-proj/R/gen_rand_zonotope.R index d4ea92eb2..3a2c8c1b2 100644 --- a/R-proj/R/gen_rand_zonotope.R +++ b/R-proj/R/gen_rand_zonotope.R @@ -1,6 +1,7 @@ #' Generator function for zonotopes #' -#' This function can be used to generate a random \eqn{d}-dimensional zonotope defined by the Minkowski sum of \eqn{m} \eqn{d}-dimensional segments. We consider \eqn{m} random directions in \eqn{R^d} and for each direction we pick a random length in \eqn{[(,\sqrt{d}]} in order to define \eqn{m} segments. +#' This function generates a random \eqn{d}-dimensional zonotope defined by the Minkowski sum of \eqn{m} \eqn{d}-dimensional segments. +#' The function considers \eqn{m} random directions in \eqn{R^d}. There are three strategies to pick the length of each segment: a) it is uniformly sampled from \eqn{[0,100]}, b) it is random from \eqn{\mathcal{N}(50,(50/3)^2)} truncated to \eqn{[0,100]}, c) it is random from \eqn{Exp(1/30)} truncated to \eqn{[0,100]}. #' #' @param dimension The dimension of the zonotope. #' @param nsegments The number of segments that generate the zonotope. diff --git a/R-proj/R/gen_simplex.R b/R-proj/R/gen_simplex.R index 9c7370605..aa2ca4a6a 100644 --- a/R-proj/R/gen_simplex.R +++ b/R-proj/R/gen_simplex.R @@ -1,6 +1,6 @@ #' Generator function for simplices #' -#' This function can be used to generate the \eqn{d}-dimensional unit simplex in H- or V-representation. +#' This function generates the \eqn{d}-dimensional unit simplex in H- or V-representation. #' #' @param dimension The dimension of the unit simplex. #' @param representation A string to declare the representation. It has to be \code{'H'} for H-representation or \code{'V'} for V-representation. diff --git a/R-proj/R/gen_skinny_cube.R b/R-proj/R/gen_skinny_cube.R index bcf3dda40..0a1df0eb0 100644 --- a/R-proj/R/gen_skinny_cube.R +++ b/R-proj/R/gen_skinny_cube.R @@ -1,6 +1,6 @@ #' Generator function for skinny hypercubes #' -#' This function can be used to generate a \eqn{d}-dimensional skinny hypercube in H-representation. +#' This function generates a \eqn{d}-dimensional skinny hypercube \eqn{[-1,1]^{d-1}\times [-100,100]}. #' #' @param dimension The dimension of the skinny hypercube. #' diff --git a/R-proj/R/round_polytope.R b/R-proj/R/round_polytope.R index 784bf90f9..246d29e54 100644 --- a/R-proj/R/round_polytope.R +++ b/R-proj/R/round_polytope.R @@ -1,12 +1,14 @@ #' Apply rounding to a convex polytope (H-polytope, V-polytope or a zonotope) #' -#' 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. +#' Given a convex H or V polytope or a zonotope as input this function brings the polytope in rounded position based on minimum volume enclosing ellipsoid of a pointset. #' #' @param P A convex polytope. It is an object from class (a) Hpolytope or (b) Vpolytope or (c) Zonotope. #' @param seed Optional. A fixed seed for the number generator. #' -#' @return A list with 2 elements: (a) a polytope of the same class as the input polytope class and (b) the element "round_value" which is the determinant of the square matrix of the linear transformation that was applied on the polytope that is given as input. +#' @return A list with 4 elements: (a) a polytope of the same class as the input polytope class and (b) the element "T" which is the matrix of the inverse linear transformation that is applied on the input polytope, (c) the element "shift" which is the opposite vector of that which has shifted the input polytope, (d) the element "round_value" which is the determinant of the square matrix of the linear transformation that is applied on the input polytope. #' +#' @references \cite{I.Z.Emiris and V. Fisikopoulos, +#' \dQuote{Practical polytope volume approximation,} \emph{ACM Trans. Math. Soft.,} 2018.}, #' @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/R/zonotope_approximation.R b/R-proj/R/zonotope_approximation.R index a66f05391..207d32e4f 100644 --- a/R-proj/R/zonotope_approximation.R +++ b/R-proj/R/zonotope_approximation.R @@ -1,24 +1,27 @@ #' A function to over-approximate a zonotope with PCA method and to evaluate the approximation by computing a ratio of fitness. #' -#' For the evaluation of the PCA method the exact volume of the approximation body is computed and the volume of the input zonotope is computed by CoolingBodies algorithm (BAN). The ratio of fitness is \eqn{R=}. +#' For the evaluation of the PCA method the exact volume of the approximation body is computed and the volume of the input zonotope is computed by CoolingBodies algorithm. The ratio of fitness is \eqn{R=vol(P) / vol(P_{red})}, where \eqn{P_{red}} is the approximate polytope. #' #' @param Z A zonotope. #' @param fit_ratio Optional. A boolean parameter to request the computation of the ratio of fitness. #' @param settings Optional. A list that declares the values of the parameters of CB algorithm as follows: #' \itemize{ -#' \item{\code{error} }{ A numeric value to set the upper bound for the approximation error. The default value is \eqn{1} for \code{'SOB'} and \eqn{0.1} otherwise.} -#' \item{\code{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{win_len} }{ The length of the sliding window for CG algorithm. The default value is \eqn{500+4dimension^2}.} -#' \item{\code{hpoly} }{ A boolean parameter to use H-polytopes in MMC of CB algorithm. The default value is \code{FALSE}.} +#' \item{\code{error} }{ A numeric value to set the upper bound for the approximation error. The default value is \eqn{0.1}.} +#' \item{\code{walk_length} }{ An integer to set the number of the steps for the random walk. The default value is \eqn{1}.} +#' \item{\code{win_len} }{ The length of the sliding window for CB algorithm. The default value is \eqn{200}.} +#' \item{\code{hpoly} }{ A boolean parameter to use H-polytopes in MMC of CB algorithm. The default value is \code{TRUE} when the order of the zonotope is \eqn{<5}, otherwise it is \code{FALSE}.} #' } #' @param seed Optional. A fixed seed for the number generator. #' #' @return A list that contains the approximation body in H-representation and the ratio of fitness #' -#' @examples +#' @references \cite{A.K. Kopetzki and B. Schurmann and M. Althoff, +#' \dQuote{Methods for Order Reduction of Zonotopes,} \emph{IEEE Conference on Decision and Control,} 2017.} +#' +#' @examples #' # over-approximate a 2-dimensional zonotope with 10 generators and compute the ratio of fitness -#' Z = gen_rand_zonotope(2,10) -#' retList = zonotope_approximation(Z = Z, fit_ratio = TRUE) +#' Z = gen_rand_zonotope(2,8) +#' retList = zonotope_approximation(Z = Z) #' #' @export zonotope_approximation <- function(Z, fit_ratio = NULL, settings = NULL, seed = NULL){ diff --git a/R-proj/man/Hpolytope.Rd b/R-proj/man/Hpolytope.Rd index 31eb2db21..c0b826747 100644 --- a/R-proj/man/Hpolytope.Rd +++ b/R-proj/man/Hpolytope.Rd @@ -3,7 +3,7 @@ \title{An \code{R} class to represent H-polytopes.} \description{ -A H-polytope is a convex polytope defined by a set of linear inequalities or equivalently a \eqn{d}-dimensional H-polytope with \eqn{m} facets is defined by a \eqn{m\times d} matrix A and a \eqn{m}-dimensional vector b, s.t.: \eqn{Ax\leq b}. +A H-polytope is a convex polytope defined by a set of linear inequalities or equivalently a \eqn{d}-dimensional H-polytope with \eqn{m} facets is defined by a \eqn{m\times d} matrix A and a \eqn{m}-dimensional vector b, s.t.: \eqn{P=\{x\ |\ Ax\leq b\} }. } \section{Fields}{ \itemize{ @@ -13,5 +13,7 @@ A H-polytope is a convex polytope defined by a set of linear inequalities or equ \item{\code{type} }{ An integer that declares the representation of the polytope. For H-representation the default value is 1.} -\item{\code{dimension} }{ An integer that declares the dimension of the polytope. -}}} +\item{\code{dimension} }{ The dimension of the polytope.} + +\item{\code{volume} }{ The volume of the polytope, if it is known.} +}} diff --git a/R-proj/man/Rcpp_Hpolytope.Rd b/R-proj/man/Rcpp_Hpolytope.Rd index 88ee372cd..88fa5b766 100644 --- a/R-proj/man/Rcpp_Hpolytope.Rd +++ b/R-proj/man/Rcpp_Hpolytope.Rd @@ -10,7 +10,7 @@ An \code{Rcpp} class to represent H-polytopes, exposed to \code{R} via modules. } \description{ -A H-polytope is a convex polytope defined by a set of linear inequalities or equivalently a \eqn{d}-dimensional H-polytope with \eqn{m} facets is defined by a \eqn{m\times d} matrix A and a \eqn{m}-dimensional vector b, s.t.: \eqn{Ax\leq b}. +A H-polytope is a convex polytope defined by a set of linear inequalities or equivalently a \eqn{d}-dimensional H-polytope with \eqn{m} facets is defined by a \eqn{m\times d} matrix A and a \eqn{m}-dimensional vector b, s.t.: \eqn{P=\{ Ax\leq b \} }. } \details{ \describe{ @@ -20,7 +20,9 @@ A H-polytope is a convex polytope defined by a set of linear inequalities or equ \item{\code{type} }{ An integer that declares the representation of the polytope. For H-representation the default value is 1.} -\item{\code{dimension} }{ An integer that declares the dimension of the polytope.} +\item{\code{dimension} }{ The dimension of the polytope.} + +\item{\code{volume} }{ The volume of the polytope, if it is known.} } } - +\keyword{internal} diff --git a/R-proj/man/Rcpp_Vpolytope.Rd b/R-proj/man/Rcpp_Vpolytope.Rd index 049d10809..95c46c099 100644 --- a/R-proj/man/Rcpp_Vpolytope.Rd +++ b/R-proj/man/Rcpp_Vpolytope.Rd @@ -18,6 +18,9 @@ A V-polytope is defined as the convex hull of \eqn{m} \eqn{d}-dimensional points \item{\code{type} }{ An integer that declares the representation of the polytope. For V-representation the default value is 2.} - \item{\code{dimension} }{ An integer that declares the dimension of the polytope.} + \item{\code{dimension} }{ The dimension of the polytope.} + + \item{\code{volume} }{ The volume of the polytope, if it is known.} } } +\keyword{internal} diff --git a/R-proj/man/Rcpp_VpolytopeIntersection.Rd b/R-proj/man/Rcpp_VpolytopeIntersection.Rd index 3ac43770b..ec308c3be 100644 --- a/R-proj/man/Rcpp_VpolytopeIntersection.Rd +++ b/R-proj/man/Rcpp_VpolytopeIntersection.Rd @@ -20,6 +20,9 @@ An intersection of two V-polytopes, \eqn{P_1}, \eqn{P_2}, is defined by the inte \item{\code{type} }{ An integer that declares the representation of the polytope. For these kinf of polytopes the default value is 4.} -\item{\code{dimension} }{ An integer that declares the dimension of the polytope.} +\item{\code{dimension} }{ The dimension of the polytope.} + +\item{\code{volume} }{ The volume of the polytope, if it is known.} } } +\keyword{internal} diff --git a/R-proj/man/Rcpp_Zonotope.Rd b/R-proj/man/Rcpp_Zonotope.Rd index 08e622d01..8c010851e 100644 --- a/R-proj/man/Rcpp_Zonotope.Rd +++ b/R-proj/man/Rcpp_Zonotope.Rd @@ -18,6 +18,9 @@ A zonotope is a convex polytope defined by the Minkowski sum of \eqn{m} \eqn{d}- \item{\code{type} }{ An integer that declares the representation of the polytope. For zonotopes the default value is 3.} -\item{\code{dimension} }{ An integer that declares the dimension of the polytope.} +\item{\code{dimension} }{ The dimension of the polytope.} + +\item{\code{volume} }{ The volume of the polytope, if it is known.} } } +\keyword{internal} diff --git a/R-proj/man/Vpolytope.Rd b/R-proj/man/Vpolytope.Rd index edcf3365b..e330af8b3 100644 --- a/R-proj/man/Vpolytope.Rd +++ b/R-proj/man/Vpolytope.Rd @@ -11,5 +11,7 @@ \item{\code{type} }{ An integer that declares the representation of the polytope. For V-representation the default value is 2.} - \item{\code{dimension} }{ An integer that declares the dimension of the polytope.} + \item{\code{dimension} }{ The dimension of the polytope.} + + \item{\code{volume} }{ The volume of the polytope, if it is known.} }} diff --git a/R-proj/man/VpolytopeIntersection.Rd b/R-proj/man/VpolytopeIntersection.Rd index 87a6801cb..f7a75c662 100644 --- a/R-proj/man/VpolytopeIntersection.Rd +++ b/R-proj/man/VpolytopeIntersection.Rd @@ -3,7 +3,7 @@ \title{An \code{R} class to represent the intersection of two V-polytopes.} \description{ - An intersection of two V-polytopes, \eqn{P_1}, \eqn{P_2}, is defined by the intersection of the two coresponding convex hulls. + An intersection of two V-polytopes, \eqn{P_1}, \eqn{P_2}, is defined by the intersection of the two corresponding convex hulls. } \section{Fields}{ \itemize{ @@ -11,7 +11,9 @@ \item{\code{V2} }{ The numerical matrix that contains the vertices of \eqn{P_2} row-wise.} -\item{\code{type} }{ An integer that declares the representation of the polytope. For these kinf of polytopes the default value is 4.} +\item{\code{type} }{ An integer that declares the representation of the polytope. For this polytope the default value is 4.} -\item{\code{dimension} }{ An integer that declares the dimension of the polytope.} +\item{\code{dimension} }{ The dimension of the polytope.} + +\item{\code{volume} }{ The volume of the polytope, if it is known.} }} diff --git a/R-proj/man/Zonotope.Rd b/R-proj/man/Zonotope.Rd index 356ebbf88..959fbed5d 100644 --- a/R-proj/man/Zonotope.Rd +++ b/R-proj/man/Zonotope.Rd @@ -11,5 +11,7 @@ A zonotope is a convex polytope defined by the Minkowski sum of \eqn{m} \eqn{d}- \item{\code{type} }{ An integer that declares the representation of the polytope. For zonotopes the default value is 3.} -\item{\code{dimension} }{ An integer that declares the dimension of the polytope.} +\item{\code{dimension} }{ The dimension of the polytope.} + +\item{\code{volume} }{ The volume of the polytope, if it is known.} }} diff --git a/R-proj/man/compute_indicators.Rd b/R-proj/man/compute_indicators.Rd index b8e01b9a4..5a79540f7 100644 --- a/R-proj/man/compute_indicators.Rd +++ b/R-proj/man/compute_indicators.Rd @@ -26,9 +26,15 @@ compute_indicators(returns, win_length = NULL, m = NULL, n = NULL, 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. } \description{ -Given a matrix that contains row-wise the assets' returns and a sliding window W, this function computes an approximation of the joint distribution (copula) between portfolios' return and volatility in each time period implied by W. -For each copula it computes an indicator: large value corresponds to a crisis period and a small value to a normal period. -The periods over which the indicator is greater than 1 for more than 60 consecutives sliding windows are warnings and for more than 100 are crisis. The sliding window is shifted by one day. +Given a matrix that contains row-wise the assets' returns and a sliding window \code{win_length}, this function computes an approximation of the joint distribution (copula, e.g. see \url{https://en.wikipedia.org/wiki/Copula_(probability_theory)}) between portfolios' return and volatility in each time period defined by \code{win_len}. +For each copula it computes an indicator: If the indicator is large it corresponds to a crisis period and if it is small it corresponds to a normal period. +In particular, the periods over which the indicator is greater than 1 for more than 60 consecutive sliding windows are warnings and for more than 100 are crisis. The sliding window is shifted by one day. +} +\examples{ +# simple example on random asset returns +asset_returns = replicate(10, rnorm(14)) +market_states_and_indicators = compute_indicators(asset_returns, 10, 10, 10000, 2, 3) + } \references{ \cite{L. Cales, A. Chalkis, I.Z. Emiris, V. Fisikopoulos, diff --git a/R-proj/man/copula.Rd b/R-proj/man/copula.Rd index 1231f54f4..4a63e1c7f 100644 --- a/R-proj/man/copula.Rd +++ b/R-proj/man/copula.Rd @@ -4,13 +4,13 @@ \alias{copula} \title{Construct a copula using uniform sampling from the unit simplex} \usage{ -copula(r1 = NULL, r2 = NULL, sigma = NULL, m = NULL, n = NULL, +copula(r1, r2 = NULL, sigma = NULL, m = NULL, n = NULL, seed = NULL) } \arguments{ -\item{r1}{A \eqn{d}-dimensional vector that describes the direction of the first family of parallel hyperplanes.} +\item{r1}{The \eqn{d}-dimensional normal vector of the first family of parallel hyperplanes.} -\item{r2}{Optional. A \eqn{d}-dimensional vector that describes the direction of the second family of parallel hyperplanes.} +\item{r2}{Optional. The \eqn{d}-dimensional normal vector of the second family of parallel hyperplanes.} \item{sigma}{Optional. The \eqn{d\times d} symmetric positive semidefine matrix that describes the family of concentric ellipsoids centered at the origin.} @@ -21,10 +21,11 @@ copula(r1 = NULL, r2 = NULL, sigma = NULL, m = NULL, n = NULL, \item{seed}{Optional. A fixed seed for the number generator.} } \value{ -A \eqn{numSlices\times numSlices} numerical matrix that corresponds to a copula. +A \eqn{m\times m} numerical matrix that corresponds to a copula. } \description{ -Given two families of parallel hyperplanes or a family of parallel hyperplanes and a family of concentric ellispoids centered at the origin intersecting the canonical simplex, this function uniformly samples from the canonical simplex and construct an approximation of the bivariate probability distribution, called copula. +Given two families of parallel hyperplanes or a family of parallel hyperplanes and a family of concentric ellispoids centered at the origin intersecting the canonical simplex, this function uniformly samples from the canonical simplex and construct an approximation of the bivariate probability distribution, called copula (see \url{https://en.wikipedia.org/wiki/Copula_(probability_theory)}). +At least two families of hyperplanes or one family of hyperplanes and one family of ellipsoids have to be given as input. } \examples{ # compute a copula for two random families of parallel hyperplanes diff --git a/R-proj/man/direct_sampling.Rd b/R-proj/man/direct_sampling.Rd index b7c9a1374..4609ec03d 100644 --- a/R-proj/man/direct_sampling.Rd +++ b/R-proj/man/direct_sampling.Rd @@ -4,7 +4,7 @@ \alias{direct_sampling} \title{Sample perfect uniformly distributed points from well known convex bodies: (a) the unit simplex, (b) the canonical simplex, (c) the boundary of a hypersphere or (d) the interior of a hypersphere.} \usage{ -direct_sampling(body = NULL, n = NULL, seed = NULL) +direct_sampling(body, n, seed = NULL) } \arguments{ \item{body}{A list to request exact uniform sampling from special well known convex bodies through the following input parameters: diff --git a/R-proj/man/exact_vol.Rd b/R-proj/man/exact_vol.Rd index b6a6a4899..593247409 100644 --- a/R-proj/man/exact_vol.Rd +++ b/R-proj/man/exact_vol.Rd @@ -19,7 +19,7 @@ For an arbitrary simplex that is given in V-representation this function compute \examples{ # compute the exact volume of a 5-dimensional zonotope defined by the Minkowski sum of 10 segments -Z = gen_rand_zonotope(5, 10) +Z = gen_rand_zonotope(2, 5) vol = exact_vol(Z) \donttest{# compute the exact volume of a 2-d arbitrary simplex @@ -32,3 +32,7 @@ vol = exact_vol(P) P = gen_cross(10,'V') vol = exact_vol(P) } +\references{ +\cite{E. Gover and N. Krikorian, +\dQuote{Determinants and the Volumes of Parallelotopes and Zonotopes,} \emph{Linear Algebra and its Applications, 433(1), 28 - 40,} 2010.} +} diff --git a/R-proj/man/file_to_polytope.Rd b/R-proj/man/file_to_polytope.Rd index dd636ccda..18f03e78b 100644 --- a/R-proj/man/file_to_polytope.Rd +++ b/R-proj/man/file_to_polytope.Rd @@ -15,7 +15,8 @@ file_to_polytope(path, zonotope) A polytope class. If the path corresponds to an ine file then the return value represents a H-polytope. If it corresponds to an ext file the return value represents a V-polytope (default choice) or a zonotope if the second argument is TRUE. } \description{ -For an ine file it generates the corresponding H-polytope. For an ext file it generates the corresponding V-polytope or zonotope. +For an ".ine" file it generates the corresponding H-polytope. For an ".ext" file it generates the corresponding V-polytope or zonotope. +For more details on those file formats see \url{https://github.com/GeomScale/volume_approximation/blob/develop/doc/cpp_interface.md#polytope-input}. } \examples{ # give the path to birk4.ine diff --git a/R-proj/man/frustum_of_simplex.Rd b/R-proj/man/frustum_of_simplex.Rd index d8b783ebd..bcb075eb4 100644 --- a/R-proj/man/frustum_of_simplex.Rd +++ b/R-proj/man/frustum_of_simplex.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/RcppExports.R \name{frustum_of_simplex} \alias{frustum_of_simplex} -\title{Compute the percentage of the volume of the unit simplex that is contained in the intersection of a half-space and the unit simplex.} +\title{Compute the percentage of the volume of the simplex that is contained in the intersection of a half-space and the simplex.} \usage{ frustum_of_simplex(a, z0) } @@ -12,10 +12,10 @@ frustum_of_simplex(a, z0) \item{z0}{The scalar that defines the half-space.} } \value{ -The percentage of the volume of the unit simplex that is contained in the intersection of a given half-space and the unit simplex. +The percentage of the volume of the simplex that is contained in the intersection of a given half-space and the simplex. } \description{ -A half-space \eqn{H} is given as a pair of a vector \eqn{a\in R^d} and a scalar \eqn{z0\in R} s.t.: \eqn{a^Tx\leq z0}. This function calls the Ali's version of the Varsi formula to compute a frustum of the unit simplex. +A half-space \eqn{H} is given as a pair of a vector \eqn{a\in R^d} and a scalar \eqn{z0\in R} s.t.: \eqn{a^Tx\leq z0}. This function calls the Ali's version of the Varsi formula to compute a frustum of the simplex. } \examples{ # compute the frustum of H: -x1+x2<=0 diff --git a/R-proj/man/gen_cross.Rd b/R-proj/man/gen_cross.Rd index 0c4df4f64..53e84c78d 100644 --- a/R-proj/man/gen_cross.Rd +++ b/R-proj/man/gen_cross.Rd @@ -15,11 +15,11 @@ gen_cross(dimension, representation) A polytope class representing a cross polytope in H- or V-representation. } \description{ -This function can be used to generate the \eqn{d}-dimensional cross polytope in H- or V-representation. +This function generates the \eqn{d}-dimensional cross polytope in H- or V-representation. } \examples{ # generate a 10-dimensional cross polytope in H-representation -P = gen_cross(10, 'H') +P = gen_cross(5, 'H') # generate a 15-dimension cross polytope in V-representation P = gen_cross(15, 'V') diff --git a/R-proj/man/gen_cube.Rd b/R-proj/man/gen_cube.Rd index b86469d2b..efb486d61 100644 --- a/R-proj/man/gen_cube.Rd +++ b/R-proj/man/gen_cube.Rd @@ -15,12 +15,12 @@ gen_cube(dimension, representation) A polytope class representing the unit \eqn{d}-dimensional hypercube in H- or V-representation. } \description{ -This function can be used to generate the \eqn{d}-dimensional unit hypercube \eqn{[-1,1]^d} in H- or V-representation. +This function generates the \eqn{d}-dimensional unit hypercube \eqn{[-1,1]^d} in H- or V-representation. } \examples{ # generate a 10-dimensional hypercube in H-representation P = gen_cube(10, 'H') # generate a 15-dimension hypercube in V-representation -P = gen_cube(15, 'V') +P = gen_cube(5, 'V') } diff --git a/R-proj/man/gen_prod_simplex.Rd b/R-proj/man/gen_prod_simplex.Rd index 3e2f599db..108a2dc1b 100644 --- a/R-proj/man/gen_prod_simplex.Rd +++ b/R-proj/man/gen_prod_simplex.Rd @@ -13,7 +13,7 @@ gen_prod_simplex(dimension) A polytope class representing the product of the two \eqn{d}-dimensional unit simplices in H-representation. } \description{ -This function can be used to generate a \eqn{2d}-dimensional polytope that is defined as the product of two \eqn{d}-dimensional unit simplices in H-representation. +This function generates a \eqn{2d}-dimensional polytope that is defined as the product of two \eqn{d}-dimensional unit simplices in H-representation. } \examples{ # generate a product of two 5-dimensional simplices. diff --git a/R-proj/man/gen_rand_hpoly.Rd b/R-proj/man/gen_rand_hpoly.Rd index b5c9adf39..e4db1c9e5 100644 --- a/R-proj/man/gen_rand_hpoly.Rd +++ b/R-proj/man/gen_rand_hpoly.Rd @@ -17,7 +17,7 @@ gen_rand_hpoly(dimension, nfacets, seed = NULL) A polytope class representing a H-polytope. } \description{ -This function can be used to generate a \eqn{d}-dimensional polytope in H-representation with \eqn{m} facets. We pick \eqn{m} random hyperplanes tangent on the \eqn{d}-dimensional unit hypersphere as facets. +This function generates a \eqn{d}-dimensional polytope in H-representation with \eqn{m} facets. We pick \eqn{m} random hyperplanes tangent on the \eqn{d}-dimensional unit hypersphere as facets. } \examples{ # generate a 10-dimensional polytope with 50 facets diff --git a/R-proj/man/gen_rand_vpoly.Rd b/R-proj/man/gen_rand_vpoly.Rd index 95dffb981..e2c3c0978 100644 --- a/R-proj/man/gen_rand_vpoly.Rd +++ b/R-proj/man/gen_rand_vpoly.Rd @@ -19,7 +19,7 @@ gen_rand_vpoly(dimension, nvertices, generator = NULL, seed = NULL) A polytope class representing a V-polytope. } \description{ -This function can be used to generate a \eqn{d}-dimensional polytope in V-representation with \eqn{m} vertices. We pick \eqn{m} random points from the boundary of the \eqn{d}-dimensional unit hypersphere as vertices. +This function generates a \eqn{d}-dimensional polytope in V-representation with \eqn{m} vertices. We pick \eqn{m} random points from the boundary of the \eqn{d}-dimensional unit hypersphere as vertices. } \examples{ # generate a 10-dimensional polytope defined as the convex hull of 25 random vertices diff --git a/R-proj/man/gen_rand_zonotope.Rd b/R-proj/man/gen_rand_zonotope.Rd index bbe3e19cc..ab971d35e 100644 --- a/R-proj/man/gen_rand_zonotope.Rd +++ b/R-proj/man/gen_rand_zonotope.Rd @@ -19,7 +19,8 @@ gen_rand_zonotope(dimension, nsegments, generator = NULL, seed = NULL) A polytope class representing a zonotope. } \description{ -This function can be used to generate a random \eqn{d}-dimensional zonotope defined by the Minkowski sum of \eqn{m} \eqn{d}-dimensional segments. We consider \eqn{m} random directions in \eqn{R^d} and for each direction we pick a random length in \eqn{[(,\sqrt{d}]} in order to define \eqn{m} segments. +This function generates a random \eqn{d}-dimensional zonotope defined by the Minkowski sum of \eqn{m} \eqn{d}-dimensional segments. +The function considers \eqn{m} random directions in \eqn{R^d}. There are three strategies to pick the length of each segment: a) it is uniformly sampled from \eqn{[0,100]}, b) it is random from \eqn{\mathcal{N}(50,(50/3)^2)} truncated to \eqn{[0,100]}, c) it is random from \eqn{Exp(1/30)} truncated to \eqn{[0,100]}. } \examples{ # generate a 10-dimensional zonotope defined by the Minkowski sum of 20 segments diff --git a/R-proj/man/gen_simplex.Rd b/R-proj/man/gen_simplex.Rd index 68b5fa086..bed04a3a6 100644 --- a/R-proj/man/gen_simplex.Rd +++ b/R-proj/man/gen_simplex.Rd @@ -15,7 +15,7 @@ gen_simplex(dimension, representation) A polytope class representing the \eqn{d}-dimensional unit simplex in H- or V-representation. } \description{ -This function can be used to generate the \eqn{d}-dimensional unit simplex in H- or V-representation. +This function generates the \eqn{d}-dimensional unit simplex in H- or V-representation. } \examples{ # generate a 10-dimensional simplex in H-representation diff --git a/R-proj/man/gen_skinny_cube.Rd b/R-proj/man/gen_skinny_cube.Rd index 9089ca1b4..916f5e20a 100644 --- a/R-proj/man/gen_skinny_cube.Rd +++ b/R-proj/man/gen_skinny_cube.Rd @@ -13,7 +13,7 @@ gen_skinny_cube(dimension) A polytope class representing the \eqn{d}-dimensional skinny hypercube in H-representation. } \description{ -This function can be used to generate a \eqn{d}-dimensional skinny hypercube in H-representation. +This function generates a \eqn{d}-dimensional skinny hypercube \eqn{[-1,1]^{d-1}\times [-100,100]}. } \examples{ # generate a 10-dimensional skinny hypercube. diff --git a/R-proj/man/inner_ball.Rd b/R-proj/man/inner_ball.Rd index 0f622d5bc..9b0c21e12 100644 --- a/R-proj/man/inner_ball.Rd +++ b/R-proj/man/inner_ball.Rd @@ -7,15 +7,14 @@ inner_ball(P) } \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) VpolytopeIntersection.} } \value{ -A \eqn{d+1}-dimensional vector that describes the inscribed ball. The first \eqn{d} coordinates corresponds to the center of the ball and the last one to the radius. +A \eqn{(d+1)}-dimensional vector that describes the inscribed ball. The first \eqn{d} coordinates corresponds to the center of the ball and the last one to the radius. } \description{ -For a H-polytope described by a \eqn{m\times d} matrix \eqn{A} and a \eqn{m}-dimensional vector \eqn{b}, s.t.: \eqn{Ax\leq b}, this function computes the largest inscribed ball (Chebychev ball) by solving the corresponding linear program. -For a V-polytope \eqn{d+1} vertices, that define a full dimensional simplex, picked at random and the largest inscribed ball of the simplex is computed. -For a zonotope \eqn{P} we compute the minimum \eqn{r} s.t.: \eqn{ r e_i \in P} for all \eqn{i=1, \dots ,d}. Then the ball centered at the origin with radius \eqn{r/ \sqrt{d}} is an inscribed ball. +For a H-polytope described by a \eqn{m\times d} matrix \eqn{A} and a \eqn{m}-dimensional vector \eqn{b}, s.t.: \eqn{P=\{x\ |\ Ax\leq b\} }, this function computes the largest inscribed ball (Chebychev ball) by solving the corresponding linear program. +For both zonotopes and V-polytopes the function computes the minimum \eqn{r} s.t.: \eqn{ r e_i \in P} for all \eqn{i=1, \dots ,d}. Then the ball centered at the origin with radius \eqn{r/ \sqrt{d}} is an inscribed ball. } \examples{ # compute the Chebychev ball of the 2d unit simplex diff --git a/R-proj/man/poly_gen.Rd b/R-proj/man/poly_gen.Rd index 30f922032..be2ed4dd0 100644 --- a/R-proj/man/poly_gen.Rd +++ b/R-proj/man/poly_gen.Rd @@ -25,8 +25,4 @@ A numerical matrix describing the requested polytope \description{ An internal Rccp function as a polytope generator } -\section{warning}{ - -Do not use this function. -} - +\keyword{internal} diff --git a/R-proj/man/rotating.Rd b/R-proj/man/rotating.Rd index 4e96293f5..4b091f815 100644 --- a/R-proj/man/rotating.Rd +++ b/R-proj/man/rotating.Rd @@ -19,8 +19,4 @@ A matrix that describes the rotated polytope \description{ An internal Rccp function for the random rotation of a convex polytope } -\section{warning}{ - -Do not use this function. -} - +\keyword{internal} diff --git a/R-proj/man/round_polytope.Rd b/R-proj/man/round_polytope.Rd index e25e90bbe..7dbce64ba 100644 --- a/R-proj/man/round_polytope.Rd +++ b/R-proj/man/round_polytope.Rd @@ -12,10 +12,10 @@ round_polytope(P, seed = NULL) \item{seed}{Optional. A fixed seed for the number generator.} } \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. +A list with 4 elements: (a) a polytope of the same class as the input polytope class and (b) the element "T" which is the matrix of the inverse linear transformation that is applied on the input polytope, (c) the element "shift" which is the opposite vector of that which has shifted the input polytope, (d) the element "round_value" which is the determinant of the square matrix of the linear transformation that is applied on the input polytope. } \description{ -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. +Given a convex H or V polytope or a zonotope as input this function brings the polytope in rounded position based on minimum volume enclosing ellipsoid of a pointset. } \examples{ # rotate a H-polytope (2d unit simplex) @@ -33,6 +33,9 @@ Z = gen_rand_zonotope(2,6) ListZono = round_polytope(Z) } \references{ +\cite{I.Z.Emiris and V. Fisikopoulos, +\dQuote{Practical polytope volume approximation,} \emph{ACM Trans. Math. Soft.,} 2018.}, + \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 ed6a6c89b..6f33c803d 100644 --- a/R-proj/man/rounding.Rd +++ b/R-proj/man/rounding.Rd @@ -12,13 +12,9 @@ rounding(P, seed = NULL) \item{seed}{Optional. A fixed seed for the number generator.} } \value{ -A numerical matrix that describes the rounded polytope and contains the round value. +A numerical matrix that describes the rounded polytope, a numerical matrix of the inverse linear transofmation that is applied on the input polytope, the numerical vector the the input polytope is shifted and the determinant of the matrix of the linear transformation that is applied on the input polytope. } \description{ Internal rcpp function for the rounding of a convex polytope } -\section{warning}{ - -Do not use this function. -} - +\keyword{internal} diff --git a/R-proj/man/sample_points.Rd b/R-proj/man/sample_points.Rd index 8cf6f5598..d4b077f62 100644 --- a/R-proj/man/sample_points.Rd +++ b/R-proj/man/sample_points.Rd @@ -2,31 +2,31 @@ % Please edit documentation in R/RcppExports.R \name{sample_points} \alias{sample_points} -\title{Sample uniformly or normally distributed points from a convex Polytope (H-polytope, V-polytope or a zonotope).} +\title{Sample uniformly or normally distributed points from a convex Polytope (H-polytope, V-polytope, zonotope or intersection of two V-polytopes).} \usage{ -sample_points(P = NULL, n = NULL, random_walk = NULL, - distribution = NULL, seed = NULL) +sample_points(P, n, random_walk = NULL, distribution = NULL, + seed = NULL) } \arguments{ -\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{P}{A convex polytope. It is an object from class (a) Hpolytope or (b) Vpolytope or (c) Zonotope or (d) VpolytopeIntersection.} \item{n}{The number of points that the function is going to sample from the convex polytope.} \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{walk_length} }{ The number of the steps per generated point for the random walk. The default value is \eqn{5} for \code{'BiW'} and \eqn{\lfloor 10 + d/10\rfloor} otherwise.} \item{\code{nburns} }{ The number of points to burn before start sampling.} -\item{\code{starting_point} }{ A \eqn{d}-dimensional numerical vector that declares a starting point in the interior of the polytope for the random walk. The default choice is the center of the Chebychev ball.} +\item{\code{starting_point} }{ A \eqn{d}-dimensional numerical vector that declares a starting point in the interior of the polytope for the random walk. The default choice is the center of the ball as that one computed by the function \code{inner_ball()}.} \item{\code{BaW_rad} }{ The radius for the ball walk.} \item{\code{L} }{ The maximum length of the billiard trajectory.} }} \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{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{mode} }{ A \eqn{d}-dimensional numerical vector that declares the mode of the Gaussian distribution. The default choice is the center of the Chebychev ball.} + \item{\code{mode} }{ A \eqn{d}-dimensional numerical vector that declares the mode of the Gaussian distribution. The default choice is the center of the as that one computed by the function \code{inner_ball()}.} }} \item{seed}{Optional. A fixed seed for the number generator.} @@ -35,11 +35,11 @@ sample_points(P = NULL, n = NULL, random_walk = NULL, A \eqn{d\times n} matrix that contains, column-wise, the sampled points from the convex polytope P. } \description{ -Sample n points with uniform or multidimensional spherical gaussian -with a mode at any point- target distribution. +Sample n points with uniform or multidimensional spherical gaussian -with a mode at any point- as the target distribution. } \examples{ -# uniform distribution from the 3d unit cube in V-representation using ball walk -P = gen_cube(3, 'V') +# uniform distribution from the 3d unit cube in H-representation using ball walk +P = gen_cube(3, 'H') 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 @@ -50,6 +50,6 @@ points = sample_points(P, n = 100, distribution = list("density" = "gaussian", " # 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")) +points = sample_points(P, n = 100, random_walk = list("walk" = "BRDHR")) } diff --git a/R-proj/man/volume.Rd b/R-proj/man/volume.Rd index 035607962..04bfd2548 100644 --- a/R-proj/man/volume.Rd +++ b/R-proj/man/volume.Rd @@ -2,21 +2,21 @@ % Please edit documentation in R/RcppExports.R \name{volume} \alias{volume} -\title{The main function for volume approximation of a convex Polytope (H-polytope, V-polytope or a zonotope)} +\title{The main function for volume approximation of a convex Polytope (H-polytope, V-polytope, zonotope or intersection of two V-polytopes)} \usage{ volume(P, settings = NULL, rounding = NULL, seed = 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) VpolytopeIntersection.} \item{settings}{Optional. A list that declares which algorithm, random walk and values of parameters to use, as follows: \itemize{ -\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{algorithm} }{ A string to set the algorithm to use: a) \code{'CB'} for CB algorithm, b) \code{'SoB'} for SOB algorithm or b) \code{'CG'} for CG algorithm. 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 SOB algorithm 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. For CB and SOB algorithms the default walk is \code{'CDHR'} for H-polytopes and \code{'BiW'} for the other representations. For CG algorithm the default walk is \code{'CDHR'} for H-polytopes and \code{'RDHR'} 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{win_len} }{ The length of the sliding window for CG algorithm. The default value is \eqn{500+4dimension^2}.} -\item{\code{hpoly} }{ A boolean parameter to use H-polytopes in MMC of CB algorithm. The default value is \code{FALSE}.} +\item{\code{win_len} }{ The length of the sliding window for CB or CG algorithm. The default value is \eqn{400+3d^2} for CB or \eqn{500+4d^2} for CG.} +\item{\code{hpoly} }{ A boolean parameter to use H-polytopes in MMC of CB algorithm when the input polytope is a zonotope. The default value is \code{TRUE} when the order of the zonotope is \eqn{<5}, otherwise it is \code{FALSE}.} }} \item{rounding}{Optional. A boolean parameter for rounding. The default value is \code{TRUE} for V-polytopes and \code{FALSE} otherwise.} @@ -27,24 +27,26 @@ volume(P, settings = NULL, rounding = NULL, seed = NULL) The approximation of the volume of a convex polytope. } \description{ -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. +For the volume approximation can be used three algorithms. Either CoolingBodies (CB) or SequenceOfBalls (SOB) or CoolingGaussian (CG). An 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{P=\{x\ |\ 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. } \examples{ -# calling SOB algorithm for a H-polytope (2d unit simplex) -P = gen_simplex(2,'H') -vol = volume(P) -# calling CG algorithm for a V-polytope (3d simplex) -P = gen_simplex(2,'V') -vol = volume(P, settings = list("algorithm" = "CG")) +# calling SOB algorithm for a H-polytope (3d unit simplex) +HP = gen_cube(3,'H') +vol = volume(HP) + +# calling CG algorithm for a V-polytope (2d simplex) +VP = gen_simplex(2,'V') +vol = volume(VP, settings = 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, settings = list("random_walk" = "RDHR", "walk_length" = 5)) +vol = volume(Z, settings = list("random_walk" = "RDHR", "walk_length" = 2)) + } \references{ \cite{I.Z.Emiris and V. Fisikopoulos, -\dQuote{Practical polytope volume approximation,} \emph{ACM Trans. Math. Soft.,} 2014.}, +\dQuote{Practical polytope volume approximation,} \emph{ACM Trans. Math. Soft.,} 2018.}, \cite{A. Chalkis and I.Z.Emiris and V. Fisikopoulos, \dQuote{Practical Volume Estimation by a New Annealing Schedule for Cooling Convex Bodies,} \emph{CoRR, abs/1905.05494,} 2019.}, diff --git a/R-proj/man/zono_approx.Rd b/R-proj/man/zono_approx.Rd index e3996dc35..9c27c369c 100644 --- a/R-proj/man/zono_approx.Rd +++ b/R-proj/man/zono_approx.Rd @@ -21,8 +21,4 @@ A List that contains a numerical matrix that describes the PCA approximation as \description{ An internal Rccp function for the over-approximation of a zonotope } -\section{warning}{ - -Do not use this function. -} - +\keyword{internal} diff --git a/R-proj/man/zonotope_approximation.Rd b/R-proj/man/zonotope_approximation.Rd index f42fdbae7..547aeb536 100644 --- a/R-proj/man/zonotope_approximation.Rd +++ b/R-proj/man/zonotope_approximation.Rd @@ -14,10 +14,10 @@ zonotope_approximation(Z, fit_ratio = NULL, settings = NULL, \item{settings}{Optional. A list that declares the values of the parameters of CB algorithm as follows: \itemize{ -\item{\code{error} }{ A numeric value to set the upper bound for the approximation error. The default value is \eqn{1} for \code{'SOB'} and \eqn{0.1} otherwise.} -\item{\code{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{win_len} }{ The length of the sliding window for CG algorithm. The default value is \eqn{500+4dimension^2}.} -\item{\code{hpoly} }{ A boolean parameter to use H-polytopes in MMC of CB algorithm. The default value is \code{FALSE}.} +\item{\code{error} }{ A numeric value to set the upper bound for the approximation error. The default value is \eqn{0.1}.} +\item{\code{walk_length} }{ An integer to set the number of the steps for the random walk. The default value is \eqn{1}.} +\item{\code{win_len} }{ The length of the sliding window for CB algorithm. The default value is \eqn{200}.} +\item{\code{hpoly} }{ A boolean parameter to use H-polytopes in MMC of CB algorithm. The default value is \code{TRUE} when the order of the zonotope is \eqn{<5}, otherwise it is \code{FALSE}.} }} \item{seed}{Optional. A fixed seed for the number generator.} @@ -26,11 +26,15 @@ zonotope_approximation(Z, fit_ratio = NULL, settings = NULL, A list that contains the approximation body in H-representation and the ratio of fitness } \description{ -For the evaluation of the PCA method the exact volume of the approximation body is computed and the volume of the input zonotope is computed by CoolingBodies algorithm (BAN). The ratio of fitness is \eqn{R=}. +For the evaluation of the PCA method the exact volume of the approximation body is computed and the volume of the input zonotope is computed by CoolingBodies algorithm. The ratio of fitness is \eqn{R=vol(P) / vol(P_{red})}, where \eqn{P_{red}} is the approximate polytope. } \examples{ # over-approximate a 2-dimensional zonotope with 10 generators and compute the ratio of fitness -Z = gen_rand_zonotope(2,10) -retList = zonotope_approximation(Z = Z, fit_ratio = TRUE) +Z = gen_rand_zonotope(2,8) +retList = zonotope_approximation(Z = Z) } +\references{ +\cite{A.K. Kopetzki and B. Schurmann and M. Althoff, +\dQuote{Methods for Order Reduction of Zonotopes,} \emph{IEEE Conference on Decision and Control,} 2017.} +} diff --git a/R-proj/src/Rproj_externals/lp_solve/lp_presolve.c b/R-proj/src/Rproj_externals/lp_solve/lp_presolve.c index 6a7dcfcba..abd7ac15d 100644 --- a/R-proj/src/Rproj_externals/lp_solve/lp_presolve.c +++ b/R-proj/src/Rproj_externals/lp_solve/lp_presolve.c @@ -32,6 +32,8 @@ initial version of column aggregation code. ------------------------------------------------------------------------- */ +//Modified by Apostolos Chalkis in May of 2020. + #include #include "commonlib.h" #include "lp_lib.h" @@ -160,7 +162,7 @@ STATIC MYBOOL presolve_fillUndo(lprec *lp, int orig_rows, int orig_cols, MYBOOL STATIC MYBOOL presolve_rebuildUndo(lprec *lp, MYBOOL isprimal) { int ik, ie, ix, j, k, *colnrDep; - LPSREAL hold, *value, *solution, *slacks; + LPSREAL hold, *value, *slacks, *solution; presolveundorec *psdata = lp->presolve_undo; MATrec *mat = NULL; @@ -168,13 +170,13 @@ STATIC MYBOOL presolve_rebuildUndo(lprec *lp, MYBOOL isprimal) if(isprimal) { if(psdata->primalundo != NULL) mat = psdata->primalundo->tracker; - solution = lp->full_solution + lp->presolve_undo->orig_rows; + //solution = lp->full_solution + lp->presolve_undo->orig_rows; // Comment out by Apostolos Chalkis slacks = lp->full_solution; } else { if(psdata->dualundo != NULL) mat = psdata->dualundo->tracker; - solution = lp->full_duals; + //solution = lp->full_duals; // Comment out by Apostolos Chalkis slacks = lp->full_duals + lp->presolve_undo->orig_rows; } if(mat == NULL) diff --git a/R-proj/src/copula.cpp b/R-proj/src/copula.cpp index eab05a2c2..065815c3e 100644 --- a/R-proj/src/copula.cpp +++ b/R-proj/src/copula.cpp @@ -19,10 +19,11 @@ //' Construct a copula using uniform sampling from the unit simplex //' -//' Given two families of parallel hyperplanes or a family of parallel hyperplanes and a family of concentric ellispoids centered at the origin intersecting the canonical simplex, this function uniformly samples from the canonical simplex and construct an approximation of the bivariate probability distribution, called copula. +//' Given two families of parallel hyperplanes or a family of parallel hyperplanes and a family of concentric ellispoids centered at the origin intersecting the canonical simplex, this function uniformly samples from the canonical simplex and construct an approximation of the bivariate probability distribution, called copula (see \url{https://en.wikipedia.org/wiki/Copula_(probability_theory)}). +//' At least two families of hyperplanes or one family of hyperplanes and one family of ellipsoids have to be given as input. //' -//' @param r1 A \eqn{d}-dimensional vector that describes the direction of the first family of parallel hyperplanes. -//' @param r2 Optional. A \eqn{d}-dimensional vector that describes the direction of the second family of parallel hyperplanes. +//' @param r1 The \eqn{d}-dimensional normal vector of the first family of parallel hyperplanes. +//' @param r2 Optional. The \eqn{d}-dimensional normal vector of the second family of parallel hyperplanes. //' @param sigma Optional. The \eqn{d\times d} symmetric positive semidefine matrix that describes the family of concentric ellipsoids centered at the origin. //' @param m The number of the slices for the copula. The default value is 100. //' @param n The number of points to sample. The default value is \eqn{5\cdot 10^5}. @@ -31,7 +32,7 @@ //' @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.} //' -//' @return A \eqn{numSlices\times numSlices} numerical matrix that corresponds to a copula. +//' @return A \eqn{m\times m} numerical matrix that corresponds to a copula. //' @examples //' # compute a copula for two random families of parallel hyperplanes //' h1 = runif(n = 10, min = 1, max = 1000) @@ -49,7 +50,7 @@ //' //' @export // [[Rcpp::export]] -Rcpp::NumericMatrix copula (Rcpp::Nullable r1 = R_NilValue, +Rcpp::NumericMatrix copula (Rcpp::Nullable r1, Rcpp::Nullable r2 = R_NilValue, Rcpp::Nullable sigma = R_NilValue, Rcpp::Nullable m = R_NilValue, @@ -79,12 +80,6 @@ Rcpp::NumericMatrix copula (Rcpp::Nullable r1 = R_NilValue, std::vector > StdCopula; unsigned int dim = Rcpp::as >(r1).size(), i, j; - if(!r1.isNotNull()) { - - throw Rcpp::exception("You have to give at least one normal of a hyperplane!"); - - } - std::vector hyp1 = Rcpp::as >(r1); if (r2.isNotNull()) { @@ -104,7 +99,7 @@ Rcpp::NumericMatrix copula (Rcpp::Nullable r1 = R_NilValue, StdCopula = hypfam_ellfam(dim, numpoints, num_slices, hyp1, Ell, seed3); } else { - throw Rcpp::exception("Wrong inputs"); + throw Rcpp::exception("You have to give as input either two families of hyperplanes or one family of hyperplanes and one family of ellipsoids!"); } diff --git a/R-proj/src/direct_sampling.cpp b/R-proj/src/direct_sampling.cpp index 8fbcbb6e1..a5a305f28 100644 --- a/R-proj/src/direct_sampling.cpp +++ b/R-proj/src/direct_sampling.cpp @@ -43,8 +43,8 @@ //' points = direct_sampling(n = 100, body = list("type" = "ball", "dimension" = 2)) //' @export // [[Rcpp::export]] -Rcpp::NumericMatrix direct_sampling(Rcpp::Nullable body = R_NilValue, - Rcpp::Nullable n = R_NilValue, +Rcpp::NumericMatrix direct_sampling(Rcpp::Nullable body, + Rcpp::Nullable n, Rcpp::Nullable seed = R_NilValue) { typedef double NT; @@ -67,7 +67,7 @@ Rcpp::NumericMatrix direct_sampling(Rcpp::Nullable body = R_NilValue RNGType rng(dim); if (seed.isNotNull()) { - unsigned seed2 = std::chrono::system_clock::now().time_since_epoch().count(); + unsigned seed2 = Rcpp::as(seed); rng.set_seed(seed2); } double seed3 = (!seed.isNotNull()) ? std::numeric_limits::signaling_NaN() : Rcpp::as(seed); @@ -75,55 +75,45 @@ Rcpp::NumericMatrix direct_sampling(Rcpp::Nullable body = R_NilValue numpoints = (!n.isNotNull()) ? 100 : Rcpp::as(n); - 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!"); - } + numpoints = Rcpp::as(n); + if (numpoints <= 0) throw Rcpp::exception("The number of samples has to be a positice integer!"); - if (body.isNotNull()) { + if (Rcpp::as(body).containsElementNamed("radius")) { - if (Rcpp::as(body).containsElementNamed("radius")) { + radius = Rcpp::as(Rcpp::as(body)["radius"]); + if (radius <= NT(0)) throw Rcpp::exception("Radius has to be a positive number!"); - radius = Rcpp::as(Rcpp::as(body)["radius"]); - if (radius <= NT(0)) throw Rcpp::exception("Radius has to be a positive number!"); + } + if (!Rcpp::as(body).containsElementNamed("type")) { - } - if (!Rcpp::as(body).containsElementNamed("type")) { + throw Rcpp::exception("The kind of body has to be given as input!"); - throw Rcpp::exception("The kind of body has to be given as input!"); + } + if (Rcpp::as(Rcpp::as(body)["type"]).compare(std::string("hypersphere"))==0) { + for (unsigned int k = 0; k < numpoints; ++k) { + randPoints.push_back(GetPointOnDsphere::apply(dim, radius, rng)); } - //RNGType rng(dim); - if (Rcpp::as(Rcpp::as(body)["type"]).compare(std::string("hypersphere"))==0) { - for (unsigned int k = 0; k < numpoints; ++k) { - randPoints.push_back(GetPointOnDsphere::apply(dim, radius, rng)); - } + } else if (Rcpp::as(Rcpp::as(body)["type"]).compare(std::string("ball"))==0) { - } else if (Rcpp::as(Rcpp::as(body)["type"]).compare(std::string("ball"))==0) { - - for (unsigned int k = 0; k < numpoints; ++k) { - randPoints.push_back(GetPointInDsphere::apply(dim, radius, rng)); - } + for (unsigned int k = 0; k < numpoints; ++k) { + randPoints.push_back(GetPointInDsphere::apply(dim, radius, rng)); + } - } else if (Rcpp::as(Rcpp::as(body)["type"]).compare(std::string("unit_simplex"))==0) { + } else if (Rcpp::as(Rcpp::as(body)["type"]).compare(std::string("unit_simplex"))==0) { - Sam_Unit(dim, numpoints, randPoints, seed3); + Sam_Unit(dim, numpoints, randPoints, seed3); - } else if (Rcpp::as(Rcpp::as(body)["type"]).compare(std::string("canonical_simplex"))==0) { + } else if (Rcpp::as(Rcpp::as(body)["type"]).compare(std::string("canonical_simplex"))==0) { - Sam_Canon_Unit(dim, numpoints, randPoints, seed3); + Sam_Canon_Unit(dim, numpoints, randPoints, seed3); - } else { + } else { - throw Rcpp::exception("Wrong input!"); + throw Rcpp::exception("Wrong input!"); - } - } else { - throw Rcpp::exception("body is null!"); } MT RetMat(dim, numpoints); diff --git a/R-proj/src/exact_vol.cpp b/R-proj/src/exact_vol.cpp index d1d3dfab8..281e9f275 100644 --- a/R-proj/src/exact_vol.cpp +++ b/R-proj/src/exact_vol.cpp @@ -31,11 +31,14 @@ FT factorial(FT n) //' //' @param P A polytope //' +//' @references \cite{E. Gover and N. Krikorian, +//' \dQuote{Determinants and the Volumes of Parallelotopes and Zonotopes,} \emph{Linear Algebra and its Applications, 433(1), 28 - 40,} 2010.} +//' //' @return The exact volume of the input polytope, for zonotopes, simplices in V-representation and polytopes with known exact volume //' @examples //' //' # compute the exact volume of a 5-dimensional zonotope defined by the Minkowski sum of 10 segments -//' Z = gen_rand_zonotope(5, 10) +//' Z = gen_rand_zonotope(2, 5) //' vol = exact_vol(Z) //' //' \donttest{# compute the exact volume of a 2-d arbitrary simplex diff --git a/R-proj/src/frustum_of_simplex.cpp b/R-proj/src/frustum_of_simplex.cpp index b9e47a86b..7841627b5 100644 --- a/R-proj/src/frustum_of_simplex.cpp +++ b/R-proj/src/frustum_of_simplex.cpp @@ -8,9 +8,9 @@ #include #include "volume/exact_vols.h" -//' Compute the percentage of the volume of the unit simplex that is contained in the intersection of a half-space and the unit simplex. +//' Compute the percentage of the volume of the simplex that is contained in the intersection of a half-space and the simplex. //' -//' A half-space \eqn{H} is given as a pair of a vector \eqn{a\in R^d} and a scalar \eqn{z0\in R} s.t.: \eqn{a^Tx\leq z0}. This function calls the Ali's version of the Varsi formula to compute a frustum of the unit simplex. +//' A half-space \eqn{H} is given as a pair of a vector \eqn{a\in R^d} and a scalar \eqn{z0\in R} s.t.: \eqn{a^Tx\leq z0}. This function calls the Ali's version of the Varsi formula to compute a frustum of the simplex. //' //' @param a A \eqn{d}-dimensional vector that defines the direction of the hyperplane. //' @param z0 The scalar that defines the half-space. @@ -21,7 +21,7 @@ //' @references \cite{Ali, Mir M., //' \dQuote{Content of the frustum of a simplex,} \emph{ Pacific J. Math. 48, no. 2, 313--322,} 1973.} //' -//' @return The percentage of the volume of the unit simplex that is contained in the intersection of a given half-space and the unit simplex. +//' @return The percentage of the volume of the simplex that is contained in the intersection of a given half-space and the simplex. //' //' @examples //' # compute the frustum of H: -x1+x2<=0 diff --git a/R-proj/src/inner_ball.cpp b/R-proj/src/inner_ball.cpp index c2c4f30a1..3c671813c 100644 --- a/R-proj/src/inner_ball.cpp +++ b/R-proj/src/inner_ball.cpp @@ -15,13 +15,12 @@ //' Compute an inscribed ball of a convex polytope //' -//' For a H-polytope described by a \eqn{m\times d} matrix \eqn{A} and a \eqn{m}-dimensional vector \eqn{b}, s.t.: \eqn{Ax\leq b}, this function computes the largest inscribed ball (Chebychev ball) by solving the corresponding linear program. -//' For a V-polytope \eqn{d+1} vertices, that define a full dimensional simplex, picked at random and the largest inscribed ball of the simplex is computed. -//' For a zonotope \eqn{P} we compute the minimum \eqn{r} s.t.: \eqn{ r e_i \in P} for all \eqn{i=1, \dots ,d}. Then the ball centered at the origin with radius \eqn{r/ \sqrt{d}} is an inscribed ball. +//' For a H-polytope described by a \eqn{m\times d} matrix \eqn{A} and a \eqn{m}-dimensional vector \eqn{b}, s.t.: \eqn{P=\{x\ |\ Ax\leq b\} }, this function computes the largest inscribed ball (Chebychev ball) by solving the corresponding linear program. +//' For both zonotopes and V-polytopes the function computes the minimum \eqn{r} s.t.: \eqn{ r e_i \in P} 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 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 or (b) Vpolytope or (c) Zonotope or (d) VpolytopeIntersection. //' -//' @return A \eqn{d+1}-dimensional vector that describes the inscribed ball. The first \eqn{d} coordinates corresponds to the center of the ball and the last one to the radius. +//' @return A \eqn{(d+1)}-dimensional vector that describes the inscribed ball. The first \eqn{d} coordinates corresponds to the center of the ball and the last one to the radius. //' //' @examples //' # compute the Chebychev ball of the 2d unit simplex diff --git a/R-proj/src/poly_gen.cpp b/R-proj/src/poly_gen.cpp index d61ccbf4a..1a8f620cf 100644 --- a/R-proj/src/poly_gen.cpp +++ b/R-proj/src/poly_gen.cpp @@ -33,8 +33,7 @@ //' @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 Optional. A fixed seed for the random polytope generator. //' -//' @section warning: -//' Do not use this function. +//' @keywords internal //' //' @return A numerical matrix describing the requested polytope // [[Rcpp::export]] diff --git a/R-proj/src/polytopes_modules.cpp b/R-proj/src/polytopes_modules.cpp index d82559a5e..e7a3dd8b8 100644 --- a/R-proj/src/polytopes_modules.cpp +++ b/R-proj/src/polytopes_modules.cpp @@ -9,67 +9,47 @@ class Hpolytope { public: Hpolytope() {} - Hpolytope(Rcpp::NumericMatrix _A, Rcpp::NumericVector _b) : A(_A), b(_b) { - dimension = _A.ncol(); - } - Hpolytope(Rcpp::NumericMatrix _A, Rcpp::NumericVector _b, double volume) : A(_A), b(_b), vol(volume) { - dimension = _A.ncol(); - } - int type = 1; - unsigned int dimension; - double vol = std::numeric_limits::signaling_NaN(); + Hpolytope(Rcpp::NumericMatrix _A, Rcpp::NumericVector _b) : A(_A), b(_b), vol(std::numeric_limits::signaling_NaN()), dimension(_A.ncol()), type(1) {} + Hpolytope(Rcpp::NumericMatrix _A, Rcpp::NumericVector _b, double volume) : A(_A), b(_b), vol(volume), dimension(_A.ncol()), type(1) {} Rcpp::NumericMatrix A; Rcpp::NumericVector b; - + double vol; + unsigned int dimension; + int type; }; class Vpolytope { public: Vpolytope() {} - Vpolytope(Rcpp::NumericMatrix _V) : V(_V) { - dimension = _V.ncol(); - } - Vpolytope(Rcpp::NumericMatrix _V, double volume) : V(_V), vol(volume) { - dimension = _V.ncol(); - } - int type = 2; - unsigned int dimension; - double vol = std::numeric_limits::signaling_NaN(); + Vpolytope(Rcpp::NumericMatrix _V) : V(_V), vol(std::numeric_limits::signaling_NaN()), dimension(_V.ncol()), type(2) {} + Vpolytope(Rcpp::NumericMatrix _V, double volume) : V(_V), vol(volume), dimension(_V.ncol()), type(2) {} Rcpp::NumericMatrix V; - + double vol; + unsigned int dimension; + int type; }; class Zonotope { public: Zonotope() {} - Zonotope(Rcpp::NumericMatrix _G) : G(_G) { - dimension = _G.ncol(); - } - Zonotope(Rcpp::NumericMatrix _G, double volume) : G(_G), vol(volume) { - dimension = _G.ncol(); - } - int type = 3; - unsigned int dimension; - double vol = std::numeric_limits::signaling_NaN(); + Zonotope(Rcpp::NumericMatrix _G) : G(_G), vol(std::numeric_limits::signaling_NaN()), dimension(_G.ncol()), type(3) {} + Zonotope(Rcpp::NumericMatrix _G, double volume) : G(_G), vol(volume), dimension(_G.ncol()), type(3) {} Rcpp::NumericMatrix G; - + double vol; + unsigned int dimension; + int type; }; class VPinterVP { public: VPinterVP() {} - VPinterVP(Rcpp::NumericMatrix _V1, Rcpp::NumericMatrix _V2) : V1(_V1), V2(_V2) { - dimension = _V1.ncol(); - } - VPinterVP(Rcpp::NumericMatrix _V1, Rcpp::NumericMatrix _V2, double volume) : V1(_V1), V2(_V2),vol(volume) { - dimension = _V1.ncol(); - } - int type = 4; - unsigned int dimension; - double vol = std::numeric_limits::signaling_NaN(); + VPinterVP(Rcpp::NumericMatrix _V1, Rcpp::NumericMatrix _V2) : V1(_V1), V2(_V2), vol(std::numeric_limits::signaling_NaN()), dimension(_V1.ncol()), type(4) {} + VPinterVP(Rcpp::NumericMatrix _V1, Rcpp::NumericMatrix _V2, double volume) : V1(_V1), V2(_V2), vol(volume), dimension(_V1.ncol()), type(4) {} Rcpp::NumericMatrix V1; Rcpp::NumericMatrix V2; - + double vol; + unsigned int dimension; + int type; }; @@ -83,9 +63,9 @@ RCPP_MODULE(yada){ //' //' @field A \eqn{m\times d} numerical matrix A //' @field b \eqn{m}-dimensional vector b - //' @field type An integer that declares the representation of the polytope. For H-representation the default value is 1. It has not be given to the constructor. - //' @field dimension An integer that declares the dimension of the polytope. It has not be given to the constructor. //' @field volume The volume of the polytope. + //' @field dimension An integer that declares the dimension of the polytope. It has not be given to the constructor. + //' @field type An integer that declares the representation of the polytope. For H-representation the default value is 1. It has not be given to the constructor. //' //' @example //' # create a 2-d unit simplex @@ -99,20 +79,21 @@ RCPP_MODULE(yada){ .constructor() .constructor() - .field( "type", &Hpolytope::type ) - .field( "dimension", &Hpolytope::dimension ) - .field( "volume", &Hpolytope::vol ) + .field( "A", &Hpolytope::A ) .field( "b", &Hpolytope::b ) - .field( "A", &Hpolytope::A ); + .field( "volume", &Hpolytope::vol ) + .field( "dimension", &Hpolytope::dimension ) + .field( "type", &Hpolytope::type ); + //' An exposed C++ class to represent a V-polytope //' //' @description A V-polytope is a convex polytope defined by the set of its vertices. //' //' @field V \eqn{m\times d} numerical matrix that contains the vertices row-wise - //' @field type An integer that declares the representation of the polytope. For H-representation the default value is 2. It has not be given to the constructor. - //' @field dimension An integer that declares the dimension of the polytope. It has not be given to the constructor. //' @field volume The volume of the polytope. + //' @field dimension An integer that declares the dimension of the polytope. It has not be given to the constructor. + //' @field type An integer that declares the representation of the polytope. For H-representation the default value is 2. It has not be given to the constructor. //' //' @example //' # Create a 2-d cube in V-representation @@ -125,19 +106,20 @@ RCPP_MODULE(yada){ .constructor() .constructor() - .field( "type", &Vpolytope::type ) - .field( "dimension", &Vpolytope::dimension ) + .field( "V", &Vpolytope::V ) .field( "volume", &Vpolytope::vol ) - .field( "V", &Vpolytope::V ); + .field( "dimension", &Vpolytope::dimension ) + .field( "type", &Vpolytope::type ); + //' An exposed C++ class to represent a zonotope //' //' @description A zonotope is a convex polytope defined by the Minkowski sum of \eqn{m} \eqn{d}-dimensional segments. //' //' @field G \eqn{m\times d} numerical matrix that contains the segments (or generators) row-wise - //' @field type An integer that declares the representation of the polytope. For H-representation the default value is 3. It has not be given to the constructor. - //' @field dimension An integer that declares the dimension of the polytope. It has not be given to the constructor. //' @field volume The volume of the polytope. + //' @field dimension An integer that declares the dimension of the polytope. It has not be given to the constructor. + //' @field type An integer that declares the representation of the polytope. For H-representation the default value is 3. It has not be given to the constructor. //' //' @example //' # Create a 2-d zonotope with 4 generators @@ -150,19 +132,21 @@ RCPP_MODULE(yada){ .constructor() .constructor() - .field( "type", &Zonotope::type ) - .field( "dimension", &Zonotope::dimension ) + .field( "G", &Zonotope::G ) .field( "volume", &Zonotope::vol ) - .field( "G", &Zonotope::G ); + .field( "dimension", &Zonotope::dimension ) + .field( "type", &Zonotope::type ); + //' An exposed C++ class to represent an intersection of two V-polytopes //' //' @description An intersection of two V-polytopes is defined by the intersection of the two coresponding convex hulls. //' - //' @field G \eqn{m\times d} numerical matrix that contains the segments (or generators) row-wise - //' @field type An integer that declares the representation of the polytope. For H-representation the default value is 4. It has not be given to the constructor. - //' @field dimension An integer that declares the dimension of the polytope. It has not be given to the constructor. + //' @field V1 \eqn{m\times d} numerical matrix that contains the vertices of the first V-polytope (row-wise) + //' @field V2 \eqn{q\times d} numerical matrix that contains the vertices of the second V-polytope (row-wise) //' @field volume The volume of the polytope. + //' @field dimension An integer that declares the dimension of the polytope. It has not be given to the constructor. + //' @field type An integer that declares the representation of the polytope. For H-representation the default value is 4. It has not be given to the constructor. //' //' @example //' # define the intwrsection of a 2-d simplex with a 2-d cross polytope @@ -176,11 +160,11 @@ RCPP_MODULE(yada){ .constructor() .constructor() - .field( "type", &VPinterVP::type ) - .field( "dimension", &VPinterVP::dimension ) - .field( "volume", &VPinterVP::vol ) .field( "V1", &VPinterVP::V1 ) - .field( "V2", &VPinterVP::V2 ); + .field( "V2", &VPinterVP::V2 ) + .field( "volume", &VPinterVP::vol ) + .field( "dimension", &VPinterVP::dimension ) + .field( "type", &VPinterVP::type ); } extern SEXP _rcpp_module_boot_yada(void); diff --git a/R-proj/src/rotating.cpp b/R-proj/src/rotating.cpp index bc86ae68c..0b688761a 100644 --- a/R-proj/src/rotating.cpp +++ b/R-proj/src/rotating.cpp @@ -25,8 +25,7 @@ //' @param T Optional. A rotation matrix. //' @param seed Optional. A fixed seed for the random linear map generator. //' -//' @section warning: -//' Do not use this function. +//' @keywords internal //' //' @return A matrix that describes the rotated polytope // [[Rcpp::export]] diff --git a/R-proj/src/rounding.cpp b/R-proj/src/rounding.cpp index 41785e503..be214d3bd 100644 --- a/R-proj/src/rounding.cpp +++ b/R-proj/src/rounding.cpp @@ -25,10 +25,9 @@ //' @param P A convex polytope (H- or V-representation or zonotope). //' @param seed Optional. A fixed seed for the number generator. //' -//' @section warning: -//' Do not use this function. +//' @keywords internal //' -//' @return A numerical matrix that describes the rounded polytope and contains the round value. +//' @return A numerical matrix that describes the rounded polytope, a numerical matrix of the inverse linear transofmation that is applied on the input polytope, the numerical vector the the input polytope is shifted and the determinant of the matrix of the linear transformation that is applied on the input polytope. // [[Rcpp::export]] Rcpp::List rounding (Rcpp::Reference P, Rcpp::Nullable seed = R_NilValue){ diff --git a/R-proj/src/sample_points.cpp b/R-proj/src/sample_points.cpp index f1bc2878d..510c74ef7 100644 --- a/R-proj/src/sample_points.cpp +++ b/R-proj/src/sample_points.cpp @@ -40,7 +40,7 @@ void sample_from_polytope(Polytope &P, RNGType &rng, PointList &randPoints, unsi uniform_sampling(randPoints, P, rng, walkL, numpoints, StartingPoint, nburns); } - } else if(rdhr){ + } else if (rdhr){ if (gaussian) { gaussian_sampling(randPoints, P, rng, walkL, numpoints, a, StartingPoint, nburns); @@ -58,7 +58,6 @@ void sample_from_polytope(Polytope &P, RNGType &rng, PointList &randPoints, unsi } } else { if (set_L) { - if (gaussian) { GaussianBallWalk WalkType(L); gaussian_sampling(randPoints, P, rng, WalkType, walkL, numpoints, a, @@ -81,33 +80,33 @@ void sample_from_polytope(Polytope &P, RNGType &rng, PointList &randPoints, unsi } -//' Sample uniformly or normally distributed points from a convex Polytope (H-polytope, V-polytope or a zonotope). +//' Sample uniformly or normally distributed points from a convex Polytope (H-polytope, V-polytope, zonotope or intersection of two V-polytopes). //' -//' Sample n points with uniform or multidimensional spherical gaussian -with a mode at any point- target distribution. +//' Sample n points with uniform or multidimensional spherical gaussian -with a mode at any point- as the target distribution. //' -//' @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 P A convex polytope. It is an object from class (a) Hpolytope or (b) Vpolytope or (c) Zonotope or (d) VpolytopeIntersection. //' @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{walk_length} }{ The number of the steps per generated point for the random walk. The default value is \eqn{5} for \code{'BiW'} and \eqn{\lfloor 10 + d/10\rfloor} otherwise.} //' \item{\code{nburns} }{ The number of points to burn before start sampling.} -//' \item{\code{starting_point} }{ A \eqn{d}-dimensional numerical vector that declares a starting point in the interior of the polytope for the random walk. The default choice is the center of the Chebychev ball.} +//' \item{\code{starting_point} }{ A \eqn{d}-dimensional numerical vector that declares a starting point in the interior of the polytope for the random walk. The default choice is the center of the ball as that one computed by the function \code{inner_ball()}.} //' \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{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{mode} }{ A \eqn{d}-dimensional numerical vector that declares the mode of the Gaussian distribution. The default choice is the center of the Chebychev ball.} +//' \item{\code{mode} }{ A \eqn{d}-dimensional numerical vector that declares the mode of the Gaussian distribution. The default choice is the center of the as that one computed by the function \code{inner_ball()}.} //' } //' @param seed Optional. A fixed seed for the number generator. //' //' @return A \eqn{d\times n} matrix that contains, column-wise, the sampled points from the convex polytope P. //' @examples -//' # uniform distribution from the 3d unit cube in V-representation using ball walk -//' P = gen_cube(3, 'V') +//' # uniform distribution from the 3d unit cube in H-representation using ball walk +//' P = gen_cube(3, 'H') //' 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 @@ -118,12 +117,12 @@ void sample_from_polytope(Polytope &P, RNGType &rng, PointList &randPoints, unsi //' //' # 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")) +//' points = sample_points(P, n = 100, random_walk = list("walk" = "BRDHR")) //' //' @export // [[Rcpp::export]] -Rcpp::NumericMatrix sample_points(Rcpp::Nullable P = R_NilValue, - Rcpp::Nullable n = R_NilValue, +Rcpp::NumericMatrix sample_points(Rcpp::Nullable P, + Rcpp::Nullable n, Rcpp::Nullable random_walk = R_NilValue, Rcpp::Nullable distribution = R_NilValue, Rcpp::Nullable seed = R_NilValue){ @@ -161,29 +160,17 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P = R_NilValue std::pair InnerBall; Point mode(dim); - numpoints = (!n.isNotNull()) ? 100 : Rcpp::as(n); - - 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 (!P.isNotNull()) { - throw Rcpp::exception("No polytope is given as input!"); - } - - + numpoints = Rcpp::as(n); + if (numpoints <= 0) throw Rcpp::exception("The number of samples has to be a positice integer!"); if (!distribution.isNotNull() || !Rcpp::as(distribution).containsElementNamed("density")) { billiard = true; } else if ( Rcpp::as(Rcpp::as(distribution)["density"]).compare(std::string("uniform")) == 0) { - billiard = true; + billiard = true; } else if ( Rcpp::as(Rcpp::as(distribution)["density"]).compare(std::string("gaussian")) == 0) { - gaussian = true; + gaussian = true; } else { throw Rcpp::exception("Wrong distribution!"); } @@ -211,15 +198,21 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P = R_NilValue if (!random_walk.isNotNull() || !Rcpp::as(random_walk).containsElementNamed("walk")) { if (gaussian) { - cdhr = true; + if (type == 1) { + cdhr = true; + } else { + rdhr = true; + } } else { billiard = true; walkL = 5; } } else if (Rcpp::as(Rcpp::as(random_walk)["walk"]).compare(std::string("CDHR")) == 0) { cdhr = true; + billiard = false; } else if (Rcpp::as(Rcpp::as(random_walk)["walk"]).compare(std::string("RDHR")) == 0) { rdhr = true; + billiard = false; } else if (Rcpp::as(Rcpp::as(random_walk)["walk"]).compare(std::string("BaW")) == 0) { if (Rcpp::as(random_walk).containsElementNamed("BaW_rad")) { L = Rcpp::as(Rcpp::as(random_walk)["BaW_rad"]); @@ -227,6 +220,7 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P = R_NilValue if (L<=0.0) throw Rcpp::exception("BaW diameter must be a postitive number!"); } ball_walk = true; + billiard = false; } 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; diff --git a/R-proj/src/volume.cpp b/R-proj/src/volume.cpp index 2bd143db2..d9d660a7a 100644 --- a/R-proj/src/volume.cpp +++ b/R-proj/src/volume.cpp @@ -80,25 +80,25 @@ double generic_volume(Polytope& P, RNGType &rng, unsigned int walk_length, NT e, return vol; } -//' The main function for volume approximation of a convex Polytope (H-polytope, V-polytope or a zonotope) +//' The main function for volume approximation of a convex Polytope (H-polytope, V-polytope, zonotope or intersection of two V-polytopes) //' -//' 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. +//' For the volume approximation can be used three algorithms. Either CoolingBodies (CB) or SequenceOfBalls (SOB) or CoolingGaussian (CG). An 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{P=\{x\ |\ 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 P A convex polytope. It is an object from class a) Hpolytope or b) Vpolytope or c) Zonotope or d) VpolytopeIntersection. //' @param settings Optional. A list that declares which algorithm, random walk and values of parameters to use, as follows: //' \itemize{ -//' \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{algorithm} }{ A string to set the algorithm to use: a) \code{'CB'} for CB algorithm, b) \code{'SoB'} for SOB algorithm or b) \code{'CG'} for CG algorithm. 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 SOB algorithm 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. For CB and SOB algorithms the default walk is \code{'CDHR'} for H-polytopes and \code{'BiW'} for the other representations. For CG algorithm the default walk is \code{'CDHR'} for H-polytopes and \code{'RDHR'} 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{win_len} }{ The length of the sliding window for CG algorithm. The default value is \eqn{500+4dimension^2}.} -//' \item{\code{hpoly} }{ A boolean parameter to use H-polytopes in MMC of CB algorithm. The default value is \code{FALSE}.} +//' \item{\code{win_len} }{ The length of the sliding window for CB or CG algorithm. The default value is \eqn{400+3d^2} for CB or \eqn{500+4d^2} for CG.} +//' \item{\code{hpoly} }{ A boolean parameter to use H-polytopes in MMC of CB algorithm when the input polytope is a zonotope. The default value is \code{TRUE} when the order of the zonotope is \eqn{<5}, otherwise it is \code{FALSE}.} //' } //' @param rounding Optional. A boolean parameter for rounding. The default value is \code{TRUE} for V-polytopes and \code{FALSE} otherwise. //' @param seed Optional. A fixed seed for the number generator. //' //' @references \cite{I.Z.Emiris and V. Fisikopoulos, -//' \dQuote{Practical polytope volume approximation,} \emph{ACM Trans. Math. Soft.,} 2014.}, +//' \dQuote{Practical polytope volume approximation,} \emph{ACM Trans. Math. Soft.,} 2018.}, //' @references \cite{A. Chalkis and I.Z.Emiris and V. Fisikopoulos, //' \dQuote{Practical Volume Estimation by a New Annealing Schedule for Cooling Convex Bodies,} \emph{CoRR, abs/1905.05494,} 2019.}, //' @references \cite{B. Cousins and S. Vempala, \dQuote{A practical volume algorithm,} \emph{Springer-Verlag Berlin Heidelberg and The Mathematical Programming Society,} 2015.} @@ -106,17 +106,19 @@ double generic_volume(Polytope& P, RNGType &rng, unsigned int walk_length, NT e, //' //' @return The approximation of the volume of a convex polytope. //' @examples -//' # calling SOB algorithm for a H-polytope (2d unit simplex) -//' P = gen_simplex(2,'H') -//' vol = volume(P) //' -//' # calling CG algorithm for a V-polytope (3d simplex) -//' P = gen_simplex(2,'V') -//' vol = volume(P, settings = list("algorithm" = "CG")) +//' # calling SOB algorithm for a H-polytope (3d unit simplex) +//' HP = gen_cube(3,'H') +//' vol = volume(HP) +//' +//' # calling CG algorithm for a V-polytope (2d simplex) +//' VP = gen_simplex(2,'V') +//' vol = volume(VP, settings = 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, settings = list("random_walk" = "RDHR", "walk_length" = 5)) +//' vol = volume(Z, settings = list("random_walk" = "RDHR", "walk_length" = 2)) +//' //' @export // [[Rcpp::export]] double volume (Rcpp::Reference P, @@ -149,7 +151,7 @@ double volume (Rcpp::Reference P, NT e; if (!Rcpp::as(settings).containsElementNamed("algorithm")) { - if (type == 2 || type == 3) { + if (type == 2 || type == 3 || type == 4) { CB = true; } else if (n <= 200) { CB = true; @@ -289,7 +291,12 @@ double volume (Rcpp::Reference P, 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 (!seed.isNotNull()) { + VPcVP.init(VP1, VP2); + } else { + unsigned seed3 = Rcpp::as(seed); + VPcVP.init(VP1, VP2, seed3); + } if (!VPcVP.is_feasible()) throw Rcpp::exception("Empty set!"); return generic_volume(VPcVP, rng, walkL, e, CG, CB, win_len, round, cdhr, rdhr, ball_walk, billiard, type); diff --git a/R-proj/src/zonotope_approximation.cpp b/R-proj/src/zonotope_approximation.cpp index ce38ee492..34c5ba255 100644 --- a/R-proj/src/zonotope_approximation.cpp +++ b/R-proj/src/zonotope_approximation.cpp @@ -24,8 +24,7 @@ //' @param settings Optional. A list that declares the values of the parameters of CB algorithm. //' @param seed Optional. A fixed seed for the number generator. //' -//' @section warning: -//' Do not use this function. +//' @keywords internal //' //' @return A List that contains a numerical matrix that describes the PCA approximation as a H-polytope and the ratio of fitness. // [[Rcpp::export]] diff --git a/R-proj/tests/testthat.R b/R-proj/tests/testthat.R index 481f3af97..051a041f8 100644 --- a/R-proj/tests/testthat.R +++ b/R-proj/tests/testthat.R @@ -1,4 +1,3 @@ library(testthat) -#library(volesti) test_check("volesti") diff --git a/R-proj/tests/testthat/test_Zvol.R b/R-proj/tests/testthat/test_Zvol.R index ffda0376c..4d1b59f4f 100644 --- a/R-proj/tests/testthat/test_Zvol.R +++ b/R-proj/tests/testthat/test_Zvol.R @@ -36,6 +36,7 @@ for (i in 1:2) { } test_that("Volume Zonotope_2_4", { + #skip_if(Sys.info()[["machine"]] %in% c("x86_32")) Z = gen_rand_zonotope(2, 4, seed = 127) res = Zruntest(Z, 'Zonotope_2_4', tol, num_of_exps, algo, 5) expect_equal(res, 1) diff --git a/R-proj/tests/testthat/test_sampling.R b/R-proj/tests/testthat/test_sampling.R index 6df5df4d9..911cf50ca 100644 --- a/R-proj/tests/testthat/test_sampling.R +++ b/R-proj/tests/testthat/test_sampling.R @@ -17,7 +17,6 @@ runsample <- function(P, name_string, dist){ } -cran_only = TRUE path = system.file('extdata', package = 'volesti') for (i in 1:2) { @@ -40,12 +39,6 @@ for (i in 1:2) { expect_equal(res, 1) }) - test_that("Sampling test", { - P = file_to_polytope(paste0(path,'/birk3.ine')) - res = runsample(P, 'H-birk3', distribution) - expect_equal(res, 1) - }) - test_that("Sampling test", { P = gen_prod_simplex(5) res = runsample(P, 'H-prod_simplex_5_5', distribution) diff --git a/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h index 0cc84bb1d..a0c7b9534 100644 --- a/include/convex_bodies/hpolytope.h +++ b/include/convex_bodies/hpolytope.h @@ -38,11 +38,7 @@ class HPolytope{ public: - HPolytope() - { - typedef typename Point::FT NT; - _inner_ball = ComputeChebychevBall(A, b); - } + HPolytope() {} std::pair InnerBall() const { diff --git a/include/convex_bodies/vpolyintersectvpoly.h b/include/convex_bodies/vpolyintersectvpoly.h index 02465d83b..076bce3e3 100644 --- a/include/convex_bodies/vpolyintersectvpoly.h +++ b/include/convex_bodies/vpolyintersectvpoly.h @@ -25,6 +25,7 @@ class IntersectionOfVpoly { typedef PointType Point; typedef typename VPolytope::MT MT; typedef typename VPolytope::VT VT; + unsigned seed; std::pair _inner_ball; NT rad; VPolytope P1; @@ -51,6 +52,13 @@ class IntersectionOfVpoly { void init(const VPolytope &P, const VPolytope &Q) { P1 = P; P2 = Q; + seed = std::chrono::system_clock::now().time_since_epoch().count(); + } + + void init(const VPolytope &P, const VPolytope &Q, unsigned &_seed) { + P1 = P; + P2 = Q; + seed = _seed; } int num_of_hyperplanes() const { @@ -116,6 +124,7 @@ class IntersectionOfVpoly { bool empty; int k = P1.get_mat().rows() + P2.get_mat().rows(); RNGType rng(k); + rng.set_seed(seed); PointInIntersection(P1.get_mat(), P2.get_mat(), GetDirection::apply(k, rng), empty); return !empty; @@ -128,6 +137,7 @@ class IntersectionOfVpoly { int k1 = V1.rows(), k2 = V2.rows(); int k = k1 + k2; RNGType rng(k); + rng.set_seed(seed); Point direction(k), p(d); std::vector vertices; typename std::vector::iterator rvert; diff --git a/include/convex_bodies/vpolytope.h b/include/convex_bodies/vpolytope.h index 7402f1eba..060607363 100644 --- a/include/convex_bodies/vpolytope.h +++ b/include/convex_bodies/vpolytope.h @@ -39,10 +39,7 @@ class VPolytope{ int *colno, *colno_mem; public: - VPolytope() - { - //_inner_ball = ComputeInnerBall(); - } + VPolytope() {} std::pair InnerBall() const { diff --git a/include/lp_oracles/misc_lp.h b/include/lp_oracles/misc_lp.h index b906af753..bbbbd92f4 100644 --- a/include/lp_oracles/misc_lp.h +++ b/include/lp_oracles/misc_lp.h @@ -1,6 +1,7 @@ #ifndef MISC_LP_H #define MISC_LP_H + #include #include #include diff --git a/include/lp_oracles/solve_lp.h b/include/lp_oracles/solve_lp.h index 6b7a4066d..d51df1250 100644 --- a/include/lp_oracles/solve_lp.h +++ b/include/lp_oracles/solve_lp.h @@ -23,6 +23,7 @@ #ifndef SOLVE_LP_H #define SOLVE_LP_H + #include #include #include diff --git a/include/lp_oracles/vpolyoracles.h b/include/lp_oracles/vpolyoracles.h index 137586832..8d31da892 100644 --- a/include/lp_oracles/vpolyoracles.h +++ b/include/lp_oracles/vpolyoracles.h @@ -1,3 +1,22 @@ +// VolEsti (volume computation and sampling library) + +// Copyright (c) 2018 Vissarion Fisikopoulos, Apostolos Chalkis + +// VolEsti is free software: you can redistribute it and/or modify it +// under the terms of the GNU Lesser General Public License as published by +// the Free Software Foundation, either version 3 of the License, or (at +// your option) any later version. +// +// VolEsti is distributed in the hope that it will be useful, but WITHOUT +// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +// for more details. +// +// See the file COPYING.LESSER for the text of the GNU Lesser General +// Public License. If you did not receive this file along with HeaDDaCHe, +// see . + + #ifndef VPOLYORACLES_H #define VPOLYORACLES_H diff --git a/include/lp_oracles/zpolyoracles.h b/include/lp_oracles/zpolyoracles.h index 79ce213cc..2a47259ce 100644 --- a/include/lp_oracles/zpolyoracles.h +++ b/include/lp_oracles/zpolyoracles.h @@ -1,3 +1,22 @@ +// VolEsti (volume computation and sampling library) + +// Copyright (c) 2018 Vissarion Fisikopoulos, Apostolos Chalkis + +// VolEsti is free software: you can redistribute it and/or modify it +// under the terms of the GNU Lesser General Public License as published by +// the Free Software Foundation, either version 3 of the License, or (at +// your option) any later version. +// +// VolEsti is distributed in the hope that it will be useful, but WITHOUT +// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +// for more details. +// +// See the file COPYING.LESSER for the text of the GNU Lesser General +// Public License. If you did not receive this file along with HeaDDaCHe, +// see . + + #ifndef ZPOLYORACLES_H #define ZPOLYORACLES_H diff --git a/include/random_walks/boundary_rdhr_walk.hpp b/include/random_walks/boundary_rdhr_walk.hpp index 185ce4fa1..7ad0ecd82 100644 --- a/include/random_walks/boundary_rdhr_walk.hpp +++ b/include/random_walks/boundary_rdhr_walk.hpp @@ -44,13 +44,15 @@ struct BRDHRWalk { for (auto j=0u; j::apply(p1.dimension(), rng); + Point v = GetDirection::apply(P.dimension(), rng); std::pair bpair = P.line_intersect(_p, v, _lamdas, _Av, _lambda); _lambda = rng.sample_urdist() * (bpair.first - bpair.second) + bpair.second; - p1 = _p + bpair.first * v; - p2 = _p + bpair.second * v; + p1 = (bpair.first * v); + p1 += _p; + p2 = (bpair.second * v); + p2 += _p; _p += (_lambda * v); } } @@ -65,7 +67,7 @@ struct BRDHRWalk _lamdas.setZero(P.num_of_hyperplanes()); _Av.setZero(P.num_of_hyperplanes()); - Point v = GetDirection::apply(p.dimension(), rng); + Point v = GetDirection::apply(P.dimension(), rng); std::pair bpair = P.line_intersect(p, v, _lamdas, _Av); _lambda = rng.sample_urdist() * (bpair.first - bpair.second) + bpair.second; _p = (_lambda * v) + p; diff --git a/include/random_walks/gaussian_ball_walk.hpp b/include/random_walks/gaussian_ball_walk.hpp index 9179d610e..de2195d94 100644 --- a/include/random_walks/gaussian_ball_walk.hpp +++ b/include/random_walks/gaussian_ball_walk.hpp @@ -46,15 +46,16 @@ struct Walk typedef typename Point::FT NT; template - static inline NT compute_delta(GenericPolytope const& P) + static inline NT compute_delta(GenericPolytope const& P, NT const& a) { - return ((P.InnerBall()).second * NT(4)) / NT(P.dimension()); + //return ((P.InnerBall()).second * NT(4)) / NT(P.dimension()); + return NT(4) * (P.InnerBall()).second / std::sqrt(std::max(NT(1), a) * NT(P.dimension())); } Walk (Polytope const& P, Point const& p, NT const& a, RandomNumberGenerator &rng) { - _delta = compute_delta(P); + _delta = compute_delta(P, a); } Walk (Polytope const& P, @@ -64,7 +65,7 @@ struct Walk parameters const& params) { _delta = params.set_delta ? params.m_L - : compute_delta(P); + : compute_delta(P, a); } template diff --git a/include/random_walks/uniform_ball_walk.hpp b/include/random_walks/uniform_ball_walk.hpp index 7db62918e..aeee06537 100644 --- a/include/random_walks/uniform_ball_walk.hpp +++ b/include/random_walks/uniform_ball_walk.hpp @@ -63,7 +63,8 @@ struct BallWalk template static inline NT compute_delta(GenericPolytope const& P) { - return ((P.InnerBall()).second * NT(4)) / NT(P.dimension()); + //return ((P.InnerBall()).second * NT(4)) / NT(P.dimension()); + return (NT(4) * (P.InnerBall()).second) / std::sqrt(NT(P.dimension())); } template diff --git a/include/random_walks/uniform_billiard_walk.hpp b/include/random_walks/uniform_billiard_walk.hpp index e704c7acd..22c4aee87 100644 --- a/include/random_walks/uniform_billiard_walk.hpp +++ b/include/random_walks/uniform_billiard_walk.hpp @@ -216,7 +216,8 @@ struct Walk template Walk(GenericPolytope const& P, Point const& p, RandomNumberGenerator &rng) { - _L = compute_diameter::template compute(P); + _Len = compute_diameter + ::template compute(P); initialize(P, p, rng); } @@ -224,7 +225,7 @@ struct Walk Walk(GenericPolytope const& P, Point const& p, RandomNumberGenerator &rng, parameters const& params) { - _L = params.set_L ? params.m_L + _Len = params.set_L ? params.m_L : compute_diameter ::template compute(P); initialize(P, p, rng); @@ -240,12 +241,12 @@ struct Walk RandomNumberGenerator &rng) { unsigned int n = P.dimension(); - NT T = rng.sample_urdist() * _L; + NT T = rng.sample_urdist() * _Len; const NT dl = 0.995; for (auto j=0u; j::apply(n, rng); Point p0 = _p; int it = 0; @@ -273,7 +274,7 @@ struct Walk inline void update_delta(NT L) { - _L = L; + _Len = L; } private : @@ -293,7 +294,7 @@ private : _p = p; _v = GetDirection::apply(n, rng); - NT T = rng.sample_urdist() * _L; + NT T = rng.sample_urdist() * _Len; Point p0 = _p; int it = 0; @@ -332,7 +333,7 @@ private : } - double _L; + NT _Len; Point _p; Point _v; NT _lambda_prev; diff --git a/include/sampling/random_point_generators.hpp b/include/sampling/random_point_generators.hpp index df99ef814..d2e5bcf37 100644 --- a/include/sampling/random_point_generators.hpp +++ b/include/sampling/random_point_generators.hpp @@ -151,7 +151,7 @@ struct BoundaryRandomPointGenerator { Walk walk(P, p, rng); Point p1(P.dimension()), p2(P.dimension()); - for (unsigned int i=0; i void uniform_sampling(PointList &randPoints, @@ -59,9 +59,9 @@ void uniform_sampling(PointList &randPoints, } template < - typename RandomNumberGenerator, typename PointList, typename Polytope, + typename RandomNumberGenerator, typename WalkTypePolicy, typename Point > @@ -97,9 +97,9 @@ void uniform_sampling(PointList &randPoints, template < typename WalkTypePolicy, - typename RandomNumberGenerator, typename PointList, typename Polytope, + typename RandomNumberGenerator, typename Point > void uniform_sampling_boundary(PointList &randPoints, @@ -124,7 +124,9 @@ void uniform_sampling_boundary(PointList &randPoints, typedef BoundaryRandomPointGenerator BoundaryRandomPointGenerator; BoundaryRandomPointGenerator::apply(P, p, nburns, walk_len, randPoints, push_back_policy, rng); + randPoints.clear(); + unsigned int n = rnum / 2; BoundaryRandomPointGenerator::apply(P, p, rnum / 2, walk_len, randPoints, push_back_policy, rng); @@ -134,9 +136,9 @@ void uniform_sampling_boundary(PointList &randPoints, template < typename WalkTypePolicy, - typename RandomNumberGenerator, typename PointList, typename Polytope, + typename RandomNumberGenerator, typename NT, typename Point > @@ -173,9 +175,9 @@ void gaussian_sampling(PointList &randPoints, template < - typename RandomNumberGenerator, typename PointList, typename Polytope, + typename RandomNumberGenerator, typename WalkTypePolicy, typename NT, typename Point diff --git a/include/sampling/simplex.hpp b/include/sampling/simplex.hpp index 2cf1c0bb1..5a42c9494 100644 --- a/include/sampling/simplex.hpp +++ b/include/sampling/simplex.hpp @@ -36,8 +36,8 @@ void Sam_Unit(unsigned int dim, 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); + unsigned rng_seed2 = seed; + rng.seed(rng_seed2); } //unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); //RNGType rng(seed); @@ -173,7 +173,6 @@ void Sam_Canon_Unit(unsigned int dim, { unsigned int j,i,x_rand,M=2147483647,pointer; // M is the largest possible integer - //std::vector x_vec; std::vector y; dim--; boost::random::uniform_int_distribution<> uidist(1,M); @@ -181,11 +180,9 @@ void Sam_Canon_Unit(unsigned int dim, 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); + unsigned rng_seed2 = seed; + rng.seed(rng_seed2); } - //unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); - //RNGType rng(seed); std::vector x_vec2; NT Ti,sum; diff --git a/include/volume/volume_cooling_balls.hpp b/include/volume/volume_cooling_balls.hpp index 2d610d86d..00c0275bc 100644 --- a/include/volume/volume_cooling_balls.hpp +++ b/include/volume/volume_cooling_balls.hpp @@ -124,95 +124,82 @@ bool get_first_ball(Polytope const& P, NT& ratio, NT const& radius_input, cooling_ball_parameters const& parameters, - RNG& rng) -{ - const unsigned max_iterarions = 10; - const unsigned tolerance = 0.00000000001; + RNG& rng) { + const unsigned max_iterarions = 20; + NT tolerance = 0.00000000001; typedef typename Polytope::PointType Point; int n = P.dimension(); int iter = 1; bool bisection_int = false; bool pass = false; bool too_few = false; - std::list randPoints; + std::list randPoints; NT rmax = parameters.rmax; NT sqrt_n = std::sqrt(NT(n)); - NT rad1 = radius_input; + NT radius1 = radius_input; - if (rmax > 0.0) - { - for (int i = 0; i < 1200; ++i) - { + if (rmax > 0.0) { + for (int i = 0; i < 1200; ++i) { randPoints.push_back(GetPointInDsphere::apply(n, rmax, rng)); } pass = check_convergence(P, randPoints, too_few, ratio, 10, true, false, parameters); - if (pass || !too_few) - { + if (pass || !too_few) { B0 = Ball(Point(n), rmax * rmax); return true; } bisection_int = true; - } else - { - rmax = 2 * sqrt_n * rad1; + } else { + rmax = 2 * sqrt_n * radius1; } - NT radius = rad1; + NT radius = radius1; - while (!bisection_int) - { + while (!bisection_int) { randPoints.clear(); too_few = false; - for (int i = 0; i < 1200; ++i) - { + for (int i = 0; i < 1200; ++i) { randPoints.push_back(GetPointInDsphere::apply(n, rmax, rng)); } if (check_convergence(P, randPoints, too_few, ratio, 10, - true, false, parameters)) - { + true, false, parameters)) { B0 = Ball(Point(n), rmax * rmax); return true; } if (too_few) break; - rad1 = rmax; + radius1 = rmax; rmax = rmax + 2 * sqrt_n * radius; } NT rad_med; - NT rad0=rad1; + NT rad0 = radius1; NT rad_m = rmax; - while (iter <= max_iterarions) - { - rad_med = 0.5*(rad1+rmax); + while (iter <= max_iterarions) { + rad_med = 0.5 * (radius1 + rmax); randPoints.clear(); too_few = false; - for (int i = 0; i < 1200; ++i) - { + for (int i = 0; i < 1200; ++i) { randPoints.push_back(GetPointInDsphere::apply(n, rad_med, rng)); } if (check_convergence(P, randPoints, too_few, ratio, 10, - true, false, parameters)) - { + true, false, parameters)) { B0 = Ball(Point(n), rad_med * rad_med); return true; } - if (too_few) - { + if (too_few) { rmax = rad_med; } else { - rad1 = rad_med; + radius1 = rad_med; } - if (rmax-rad1 < tolerance) - { - rad1 = rad0; + if (rmax - radius1 < tolerance) { + radius1 = rad0; rmax = rad_m; iter++; } @@ -228,8 +215,8 @@ bool get_next_zonotopeball(std::vector& BallSet, std::vector& ratios, cooling_ball_parameters const& parameters) { - const unsigned max_iterarions = 10; - const unsigned tolerance = 0.00000000001; + const unsigned max_iterarions = 20; + NT tolerance = 0.00000000001; int n = (*randPoints.begin()).dimension(); int iter = 1; bool too_few; @@ -703,8 +690,9 @@ NT estimate_ratio_interval(PolyBall1 const& Pb1, template < typename WalkTypePolicy, - typename RandomNumberGenerator, - typename Polytope + typename Polytope, + typename RandomNumberGenerator + > double volume_cooling_balls(Polytope const& Pin, RandomNumberGenerator &rng, diff --git a/include/volume/volume_cooling_gaussians.hpp b/include/volume/volume_cooling_gaussians.hpp index 16f6fba1c..9b12b4d1c 100644 --- a/include/volume/volume_cooling_gaussians.hpp +++ b/include/volume/volume_cooling_gaussians.hpp @@ -305,8 +305,9 @@ struct gaussian_annealing_parameters template < typename WalkTypePolicy, - typename RandomNumberGenerator, - typename Polytope + typename Polytope, + typename RandomNumberGenerator + > double volume_cooling_gaussians(Polytope const& Pin, RandomNumberGenerator& rng, @@ -450,7 +451,8 @@ double volume_cooling_gaussians(Polytope const& Pin, done=true; } - index = index==W ? 0 : index%W + 1; + index = index%W + 1; + if (index == W) index = 0; } #ifdef VOLESTI_DEBUG std::cout << "ratio " << i << " = " << (*fnIt) / (*itsIt) diff --git a/include/volume/volume_cooling_hpoly.hpp b/include/volume/volume_cooling_hpoly.hpp index 1dbc02de2..dd154c67e 100644 --- a/include/volume/volume_cooling_hpoly.hpp +++ b/include/volume/volume_cooling_hpoly.hpp @@ -11,8 +11,6 @@ #include "volume/volume_cooling_gaussians.hpp" #include "sampling/random_point_generators.hpp" -#define MAX_ITER 20 -#define TOL 0.00000000001 template < @@ -36,6 +34,9 @@ bool get_first_poly(Zonotope &P, PushBackWalkPolicy push_back_policy; + const unsigned max_iterarions = 20; + NT tolerance = 0.00000000001; + MT G = P.get_mat().transpose(), A = HP.get_mat(); b_max = (A*G).cwiseAbs().rowwise().sum(); VT b_min = HP.get_vec(); @@ -49,7 +50,7 @@ bool get_first_poly(Zonotope &P, NT l=0.0, u=1.0, med; VT b_med(m); - while(iter <= MAX_ITER) { + while(iter <= max_iterarions) { q=Point(n); med = (u + l) * 0.5; @@ -79,7 +80,7 @@ bool get_first_poly(Zonotope &P, HP.set_vec(b_med); return true; } - if(u-l < TOL) { + if(u-l < tolerance) { u=1.0; l=0.0; iter++; @@ -101,12 +102,15 @@ bool get_next_zonoball(std::vector &HPolySet, typedef typename Zonotope::PointType Point; + const unsigned max_iterarions = 20; + NT tolerance = 0.00000000001; + int n = HP2.dimension(), iter = 1; bool too_few; VT b_med(b_max.size()); NT ratio, med, u = 1.0, l = 0.0; - while (iter <= MAX_ITER) { + while (iter <= max_iterarions) { med = (u + l) * 0.5; b_med = b_min + (b_max-b_min) * med; HP2.set_vec(b_med); @@ -124,7 +128,7 @@ bool get_next_zonoball(std::vector &HPolySet, } else { u = med; } - if(u-l < TOL) { + if(u-l < tolerance) { u=1.0; l=0.0; iter++; @@ -322,7 +326,7 @@ double volume_cooling_hpoly (Zonotope const& Pin, 10 + 10 * n, rng); //TODO: rounding to HP2 //NT vol = res.second * volume_cooling_balls(HP2, Her, 1); - NT vol = res.second * volume_cooling_gaussians(HP2, Her/2.0, 1); + NT vol = res.second * volume_cooling_gaussians(HP2, rng, Her/2.0, 1); if (!parameters.window2) { vol *= estimate_ratio_interval(HP, P, ratio, er0, parameters.win_len, 1200, prob, 10+10*n, rng); diff --git a/include/volume/volume_sequence_of_balls.hpp b/include/volume/volume_sequence_of_balls.hpp index da3058d8d..d3c0050f9 100644 --- a/include/volume/volume_sequence_of_balls.hpp +++ b/include/volume/volume_sequence_of_balls.hpp @@ -43,8 +43,9 @@ template < typename WalkTypePolicy, - typename RandomNumberGenerator, - typename Polytope + typename Polytope, + typename RandomNumberGenerator + > double volume_sequence_of_balls(Polytope const& Pin, RandomNumberGenerator &rng, diff --git a/test/volume_cb_hpolytope.cpp b/test/volume_cb_hpolytope.cpp index c0eeb0d33..d36de032a 100644 --- a/test/volume_cb_hpolytope.cpp +++ b/test/volume_cb_hpolytope.cpp @@ -81,11 +81,11 @@ void call_test_cube(){ std::cout << "--- Testing volume of H-cube10" << std::endl; P = gen_cube(10, false); - test_volume(P, 1151.18, 1163.36, 1119.15, 1100.73, 1024); + test_volume(P, 1118.63, 1163.36, 1119.15, 1100.73, 1024); std::cout << "--- Testing volume of H-cube20" << std::endl; P = gen_cube(20, false); - test_volume(P, 968319, 1051230, 1006470, 1007020, 1048576); + test_volume(P, 965744, 1051230, 1006470, 1007020, 1048576); } template @@ -116,7 +116,7 @@ void call_test_cross(){ std::cout << "--- Testing volume of H-cross10" << std::endl; Hpolytope P = gen_cross(10, false); test_volume(P, - 0.000296924, + 0.000291034, 0.000281135, 0.000293788, 0.000286311, @@ -139,7 +139,7 @@ void call_test_birk() inp.open("../R-proj/inst/extdata/birk3.ine",std::ifstream::in); read_pointset(inp,Pin); P.init(Pin); - test_volume(P, 0.111157, 0.125548, 0.113241, 0.116259, 0.125); + test_volume(P, 0.114343, 0.125548, 0.113241, 0.116259, 0.125); std::cout << "--- Testing volume of H-birk4" << std::endl; std::ifstream inp2; @@ -147,7 +147,7 @@ void call_test_birk() inp2.open("../R-proj/inst/extdata/birk4.ine",std::ifstream::in); read_pointset(inp2,Pin2); P.init(Pin2); - test_volume(P, 0.00074316, 0.00109593, 0.000881856, 0.000839499, + test_volume(P, 0.00106935, 0.00109593, 0.000881856, 0.000839499, 0.000970018); std::cout << "--- Testing volume of H-birk5" << std::endl; @@ -157,7 +157,7 @@ void call_test_birk() read_pointset(inp3,Pin3); P.init(Pin3); test_volume(P, - 4.1826 * std::pow(10,-9), + 9.47562 * std::pow(10,-8), 2.12236 * std::pow(10,-7), 1.87499 * std::pow(10,-7), 1.93315 * std::pow(10,-7), @@ -170,7 +170,7 @@ void call_test_birk() read_pointset(inp4,Pin4); P.init(Pin4); test_volume(P, - 3.43196 * std::pow(10,-19), + 1.95177 * std::pow(10,-14), 6.60745 * std::pow(10,-13), 5.99551 * std::pow(10,-13), 9.81049 * std::pow(10,-13), @@ -188,7 +188,7 @@ void call_test_prod_simplex() { std::cout << "--- Testing volume of H-prod_simplex5" << std::endl; P = gen_prod_simplex(5); test_volume(P, - 8.00834 * std::pow(10,-5), + 6.40072 * std::pow(10,-5), 6.69062 * std::pow(10,-5), 7.44088 * std::pow(10,-5), 6.31986 * std::pow(10,-5), @@ -197,7 +197,7 @@ void call_test_prod_simplex() { std::cout << "--- Testing volume of H-prod_simplex10" << std::endl; P = gen_prod_simplex(10); test_volume(P, - 3.77706 * std::pow(10,-14), + 6.83631 * std::pow(10,-14), 8.19581 * std::pow(10,-14), 7.42207 * std::pow(10,-14), 8.1113 * std::pow(10,-14), @@ -206,7 +206,7 @@ void call_test_prod_simplex() { std::cout << "--- Testing volume of H-prod_simplex15" << std::endl; P = gen_prod_simplex(15); test_volume(P, - 1.97215 * std::pow(10,-26), + 6.25978 * std::pow(10,-25), 9.33162 * std::pow(10,-25), 6.01102 * std::pow(10,-25), 6.45706 * std::pow(10,-25), @@ -224,7 +224,7 @@ void call_test_simplex() { std::cout << "--- Testing volume of H-simplex10" << std::endl; P = gen_simplex(10, false); test_volume(P, - 2.82928 * std::pow(10,-7), + 3.22432 * std::pow(10,-7), 2.90617 * std::pow(10,-7), 2.93392 * std::pow(10,-7), 3.03629 * std::pow(10,-7), @@ -233,7 +233,7 @@ void call_test_simplex() { std::cout << "--- Testing volume of H-simplex20" << std::endl; P = gen_simplex(20, false); test_volume(P, - 1.61574 * std::pow(10,-19), + 4.03788 * std::pow(10,-19), 4.14182 * std::pow(10,-19), 3.8545 * std::pow(10,-19), 4.28227 * std::pow(10,-19), @@ -242,7 +242,7 @@ void call_test_simplex() { std::cout << "--- Testing volume of H-simplex30" << std::endl; P = gen_simplex(30, false); test_volume(P, - 1.00967 * std::pow(10,-33), + 2.5776 * std::pow(10,-33), 3.5157 * std::pow(10,-33), 3.53407 * std::pow(10,-33), 3.59591 * std::pow(10,-33), diff --git a/test/volume_cb_vpoly_intersection_vpoly.cpp b/test/volume_cb_vpoly_intersection_vpoly.cpp index 041e6cb30..f2cd65ce5 100644 --- a/test/volume_cb_vpoly_intersection_vpoly.cpp +++ b/test/volume_cb_vpoly_intersection_vpoly.cpp @@ -62,26 +62,27 @@ void test_volume(Polytope &P1, Polytope &P2, std::cout << "Number type: " << typeid(NT).name() << std::endl; typedef BoostRandomNumberGenerator RNGType; + unsigned seed = 105; //TODO: low accuracy in high dimensions - P.init(P1, P2); + P.init(P1, P2, seed); NT volume = volume_cooling_balls(P, e/2.0, walk_len); test_values(volume, expectedBall, exact); P1.init(P.dimension(), P1.get_mat(), P1.get_vec()); P2.init(P.dimension(), P2.get_mat(), P2.get_vec()); - P.init(P1, P2); + P.init(P1, P2, seed); volume = volume_cooling_balls(P, e, walk_len); test_values(volume, expectedCDHR, exact); P1.init(P.dimension(), P1.get_mat(), P1.get_vec()); P2.init(P.dimension(), P2.get_mat(), P2.get_vec()); - P.init(P1, P2); + P.init(P1, P2, seed); volume = volume_cooling_balls(P, e, walk_len); test_values(volume, expectedRDHR, exact); P1.init(P.dimension(), P1.get_mat(), P1.get_vec()); P2.init(P.dimension(), P2.get_mat(), P2.get_vec()); - P.init(P1, P2); + P.init(P1, P2, seed); volume = volume_cooling_balls(P, e, walk_len); test_values(volume, expectedBilliard, exact); } diff --git a/test/volume_cb_zonotopes.cpp b/test/volume_cb_zonotopes.cpp index 202939c32..2ac57c406 100644 --- a/test/volume_cb_zonotopes.cpp +++ b/test/volume_cb_zonotopes.cpp @@ -129,7 +129,7 @@ void call_test_uniform_generator(){ exact_vol = exact_zonotope_vol(P); test_volume_hpoly(P, 0, - 6.35348 * std::pow(10,20), + 7.1588 * std::pow(10,20), 6.46196 * std::pow(10,20), 5.98586 * std::pow(10,20), exact_vol); diff --git a/test/volume_cg_hpolytope.cpp b/test/volume_cg_hpolytope.cpp index 3cd5e73f0..70ccabe95 100644 --- a/test/volume_cg_hpolytope.cpp +++ b/test/volume_cg_hpolytope.cpp @@ -77,11 +77,11 @@ void call_test_cube(){ std::cout << "--- Testing volume of H-cube10" << std::endl; P = gen_cube(10, false); - test_volume(P, 1102.47, 1119.38, 1037.34, 1024); + test_volume(P, 1079.56, 1110.92, 1113.93, 1024); std::cout << "--- Testing volume of H-cube20" << std::endl; P = gen_cube(20, false); - test_volume(P, 1104980, 1051120, 989794, 1048576); + test_volume(P, 1.1025e+06, 1.05174e+06, 995224, 1048576); } template @@ -112,9 +112,9 @@ void call_test_cross(){ std::cout << "--- Testing volume of H-cross10" << std::endl; Hpolytope P = gen_cross(10, false); test_volume(P, - 0.000280513, - 0.000274237, - 0.000296688, + 0.000292199, + 0.000274014, + 0.000294463, 0.0002821869); } @@ -133,7 +133,7 @@ void call_test_birk() { inp.open("../R-proj/inst/extdata/birk3.ine",std::ifstream::in); read_pointset(inp,Pin); P.init(Pin); - test_volume(P, 0.101546, 0.121466, 0.117802, 0.125); + test_volume(P, 0.116678, 0.122104, 0.11326, 0.125); std::cout << "--- Testing volume of H-birk4" << std::endl; std::ifstream inp2; @@ -142,9 +142,9 @@ void call_test_birk() { read_pointset(inp2,Pin2); P.init(Pin2); test_volume(P, - 0.000942906, - 0.00114121, - 0.00104026, + 0.000450761, + 0.00108943, + 0.00103573, 0.000970018); std::cout << "--- Testing volume of H-birk5" << std::endl; @@ -154,9 +154,9 @@ void call_test_birk() { read_pointset(inp3,Pin3); P.init(Pin3); test_volume(P, - 1.08184 * std::pow(10,-7), - 2.31411 * std::pow(10,-7), - 2.39421 * std::pow(10,-7), + 2.97522 * std::pow(10,-8), + 2.25982 * std::pow(10,-7), + 2.24768 * std::pow(10,-7), 2.25 * std::pow(10,-7)); std::cout << "--- Testing volume of H-birk6" << std::endl; @@ -166,9 +166,9 @@ void call_test_birk() { read_pointset(inp4,Pin4); P.init(Pin4); test_volume(P, - 5.35067 * std::pow(10,-17), - 8.44699 * std::pow(10,-13), - 9.58445 * std::pow(10,-13), + 3.66375 * std::pow(10,-19), + 9.85929 * std::pow(10,-13), + 1.05038 * std::pow(10,-12), 9.455459196 * std::pow(10,-13)); } @@ -183,25 +183,25 @@ void call_test_prod_simplex() { std::cout << "--- Testing volume of H-prod_simplex5" << std::endl; P = gen_prod_simplex(5); test_volume(P, - 8.61492 * std::pow(10,-5), - 7.36265 * std::pow(10,-5), - 7.07551 * std::pow(10,-5), + 6.3448 * std::pow(10,-5), + 6.94695 * std::pow(10,-5), + 6.2862 * std::pow(10,-5), std::pow(1.0 / factorial(5.0), 2)); std::cout << "--- Testing volume of H-prod_simplex10" << std::endl; P = gen_prod_simplex(10); test_volume(P, - 8.9345 * std::pow(10,-14), - 8.49473 * std::pow(10,-14), - 7.36811 * std::pow(10,-14), + 1.36206 * std::pow(10,-14), + 8.48116 * std::pow(10,-14), + 7.37191 * std::pow(10,-14), std::pow(1.0 / factorial(10.0), 2)); std::cout << "--- Testing volume of H-prod_simplex15" << std::endl; P = gen_prod_simplex(15); test_volume(P, - 6.37869 * std::pow(10,-25), - 5.55858 * std::pow(10,-25), - 5.93187 * std::pow(10,-25), + 1.93763 * std::pow(10,-26), + 5.4624 * std::pow(10,-25), + 4.99596 * std::pow(10,-25), std::pow(1.0 / factorial(15.0), 2)); } @@ -216,25 +216,25 @@ void call_test_simplex() { std::cout << "--- Testing volume of H-simplex10" << std::endl; P = gen_simplex(10, false); test_volume(P, - 3.0842 * std::pow(10,-7), - 2.59229 * std::pow(10,-7), - 2.68156 * std::pow(10,-7), + 3.14369 * std::pow(10,-7), + 2.70598 * std::pow(10,-7), + 2.87717 * std::pow(10,-7), 1.0 / factorial(10.0)); std::cout << "--- Testing volume of H-simplex20" << std::endl; P = gen_simplex(20, false); test_volume(P, - 3.41016 * std::pow(10,-19), - 4.16827 * std::pow(10,-19), - 4.43788 * std::pow(10,-19), + 1.12891 * std::pow(10,-23), + 4.16845 * std::pow(10,-19), + 3.93817 * std::pow(10,-19), 1.0 / factorial(20.0)); std::cout << "--- Testing volume of H-simplex30" << std::endl; P = gen_simplex(30, false); test_volume(P, - 3.30097 * std::pow(10,-33), - 4.15532 * std::pow(10,-33), - 3.54527 * std::pow(10,-33), + 7.06547 * std::pow(10,-41), + 4.02288 * std::pow(10,-33), + 3.41361 * std::pow(10,-33), 1.0 / factorial(30.0)); } diff --git a/test/volume_cg_vpolytope.cpp b/test/volume_cg_vpolytope.cpp index 50a57be52..826fd45e6 100644 --- a/test/volume_cg_vpolytope.cpp +++ b/test/volume_cg_vpolytope.cpp @@ -123,17 +123,17 @@ void call_test_cross(){ std::cout << "--- Testing volume of V-cross5" << std::endl; P = gen_cross(5, true); test_volume(P, - 0.274712, - 0.27768, - 0.251894, + 0.274801, + 0.277746, + 0.251418, 0.266666667); std::cout << "--- Testing volume of V-cross10" << std::endl; P = gen_cross(10, true); test_volume(P, - 0.000288736, - 0.000306377, - 0.000296662, + 0.000309838, + 0.000311191, + 0.000299492, 0.0002821869); // both slow and inaccurate for CG // std::cout << "--- Testing volume of V-cross20" << std::endl; @@ -157,9 +157,9 @@ void call_test_simplex() { std::cout << "--- Testing volume of V-simplex5" << std::endl; P = gen_simplex(5, true); test_volume(P, - 0.00845454, - 0.00855895, - 0.00837661, + 0.00694196, + 0.00885402, + 0.00842199, 1.0 / factorial(5.0)); // too slow for CG /* diff --git a/test/volume_sob_hpolytope.cpp b/test/volume_sob_hpolytope.cpp index 01b2730d0..84256bb0f 100644 --- a/test/volume_sob_hpolytope.cpp +++ b/test/volume_sob_hpolytope.cpp @@ -83,11 +83,11 @@ void call_test_cube(){ std::cout << "--- Testing volume of H-cube10" << std::endl; P = gen_cube(10, false); - test_volume(P, 1096.5089688155, 1049.22, 1055.73, 1055.73, 1024); + test_volume(P, 1014.69, 1049.22, 1055.73, 1055.73, 1024); std::cout << "--- Testing volume of H-cube20" << std::endl; P = gen_cube(20, false); - test_volume(P, 967353, 1056180, 1058830, 1058830, 1048576); + test_volume(P, 1029780, 1056180, 1058830, 1058830, 1048576); } template @@ -123,7 +123,7 @@ void call_test_cross(){ std::cout << "--- Testing volume of H-cross10" << std::endl; Hpolytope P = gen_cross(10, false); test_volume(P, - 0.000280621, + 0.000283788, 0.000280815, 0.000296745, 0.000296745, @@ -146,7 +146,7 @@ void call_test_birk() { read_pointset(inp,Pin); P.init(Pin); test_volume(P, - 0.118885, + 0.130806, 0.126776, 0.122177, 0.122177, @@ -159,7 +159,7 @@ void call_test_birk() { read_pointset(inp2,Pin2); P.init(Pin2); test_volume(P, - 0.00122254, + 0.00112925, 0.000898527, 0.000945447, 0.000945447, @@ -172,7 +172,7 @@ void call_test_birk() { read_pointset(inp3,Pin3); P.init(Pin3); test_volume(P, - 2.19189 * std::pow(10,-8), + 1.8241 * std::pow(10,-7), 2.07943 * std::pow(10,-7), 2.28319 * std::pow(10,-7), 2.28319 * std::pow(10,-7), @@ -185,7 +185,7 @@ void call_test_birk() { read_pointset(inp4,Pin4); P.init(Pin4); test_volume(P, - 5.72936 * std::pow(10,-18), + 5.27883 * std::pow(10,-13), 9.48912 * std::pow(10,-13), 7.05416 * std::pow(10,-13), 7.05416 * std::pow(10,-13), @@ -203,7 +203,7 @@ void call_test_prod_simplex() { std::cout << "--- Testing volume of H-prod_simplex5" << std::endl; P = gen_prod_simplex(5); test_volume(P, - 6.21239 * std::pow(10,-5), + 7.35223 * std::pow(10,-5), 6.86576 * std::pow(10,-5), 7.43136 * std::pow(10,-5), 7.43136 * std::pow(10,-5), @@ -212,7 +212,7 @@ void call_test_prod_simplex() { std::cout << "--- Testing volume of H-prod_simplex10" << std::endl; P = gen_prod_simplex(10); test_volume(P, - 5.92854 * std::pow(10,-14), + 7.38925 * std::pow(10,-14), 8.01351 * std::pow(10,-14), 8.27387 * std::pow(10,-14), 8.27387 * std::pow(10,-14), @@ -221,7 +221,7 @@ void call_test_prod_simplex() { std::cout << "--- Testing volume of H-prod_simplex15" << std::endl; P = gen_prod_simplex(15); test_volume(P, - 6.06129 * std::pow(10,-25), + 5.61238 * std::pow(10,-25), 5.87558 * std::pow(10,-25), 5.48179 * std::pow(10,-25), 5.48179 * std::pow(10,-25), @@ -239,7 +239,7 @@ void call_test_simplex() { std::cout << "--- Testing volume of H-simplex10" << std::endl; P = gen_simplex(10, false); test_volume(P, - 2.99056 * std::pow(10,-7), + 2.98074 * std::pow(10,-7), 2.52756 * std::pow(10,-7), 2.89366 * std::pow(10,-7), 2.89366 * std::pow(10,-7), @@ -248,7 +248,7 @@ void call_test_simplex() { std::cout << "--- Testing volume of H-simplex20" << std::endl; P = gen_simplex(20, false); test_volume(P, - 3.181 * std::pow(10,-19), + 4.64611 * std::pow(10,-19), 4.4317 * std::pow(10,-19), 4.16737 * std::pow(10,-19), 4.16737 * std::pow(10,-19), @@ -257,7 +257,7 @@ void call_test_simplex() { std::cout << "--- Testing volume of H-simplex30" << std::endl; P = gen_simplex(30, false); test_volume(P, - 2.4831 * std::pow(10,-33), + 3.65853 * std::pow(10,-33), 3.86474 * std::pow(10,-33), 4.04136 * std::pow(10,-33), 4.04136 * std::pow(10,-33), diff --git a/test/volume_sob_vpolytope.cpp b/test/volume_sob_vpolytope.cpp index 99ffb91f0..2e0276d67 100644 --- a/test/volume_sob_vpolytope.cpp +++ b/test/volume_sob_vpolytope.cpp @@ -99,11 +99,11 @@ void call_test_cross(){ std::cout << "--- Testing volume of V-cross5" << std::endl; P = gen_cross(5, true); - test_volume(P, 0.271736, 0.266666667); + test_volume(P, 0.265422, 0.266666667); std::cout << "--- Testing volume of V-cross10" << std::endl; P = gen_cross(10, true); - test_volume(P, 0.000280621, 0.0002821869); + test_volume(P, 0.000283788, 0.0002821869); } template @@ -117,7 +117,7 @@ void call_test_simplex() { std::cout << "--- Testing volume of V-simplex5" << std::endl; P = gen_simplex(5, true); - test_volume(P, 0.0100291, 1.0 / factorial(10.0)); + test_volume(P, 0.00827958, 1.0 / factorial(10.0)); // too slow for SoB // std::cout << "--- Testing volume of V-simplex10" << std::endl;