R/RcppExports.R

In volesti: Volume Approximation and Sampling of Convex Polytopes

# Generated by using Rcpp::compileAttributes() -> do not edit by hand
# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393

#' 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 (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 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}.
#' @param seed Optional. A fixed seed for the number generator.
#'
#' @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{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)
#' h1 = h1 / 1000
#' h2=runif(n = 10, min = 1, max = 1000)
#' h2 = h2 / 1000
#' cop = copula(r1 = h1, r2 = h2, m = 10, n = 100000)
#'
#' # compute a copula for a family of parallel hyperplanes and a family of conentric ellipsoids
#' h = runif(n = 10, min = 1, max = 1000)
#' h = h / 1000
#' E = replicate(10, rnorm(20))
#' E = cov(E)
#' cop = copula(r1 = h, sigma = E, m = 10, n = 100000)
#'
#' @export
copula <- function(r1, r2 = NULL, sigma = NULL, m = NULL, n = NULL, seed = NULL) {
    .Call(`_volesti_copula`, r1, r2, sigma, m, n, seed)
}

#' 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.
#'
#' The \eqn{d}-dimensional unit simplex is the set of points \eqn{\vec{x}\in \R^d}, s.t.: \eqn{\sum_i x_i\leq 1}, \eqn{x_i\geq 0}. The \eqn{d}-dimensional canonical simplex is the set of points \eqn{\vec{x}\in \R^d}, s.t.: \eqn{\sum_i x_i = 1}, \eqn{x_i\geq 0}.
#'
#' @param body A list to request exact uniform sampling from special well known convex bodies through the following input parameters:
#' \itemize{
#' \item{\code{type} }{ A string that declares the type of the body for the exact sampling: a) \code{'unit_simplex'} for the unit simplex, b) \code{'canonical_simplex'} for the canonical simplex, c) \code{'hypersphere'} for the boundary of a hypersphere centered at the origin, d) \code{'ball'} for the interior of a hypersphere centered at the origin.}
#' \item{\code{dimension} }{ An integer that declares the dimension when exact sampling is enabled for a simplex or a hypersphere.}
#' \item{\code{radius} }{ The radius of the \eqn{d}-dimensional hypersphere. The default value is \eqn{1}.}
#' \item{\code{seed} }{ A fixed seed for the number generator.}
#' }
#' @param n The number of points that the function is going to sample.
#'
#' @references \cite{R.Y. Rubinstein and B. Melamed,
#' \dQuote{Modern simulation and modeling} \emph{ Wiley Series in Probability and Statistics,} 1998.}
#' @references \cite{A Smith, Noah and W Tromble, Roy,
#' \dQuote{Sampling Uniformly from the Unit Simplex,} \emph{ Center for Language and Speech Processing Johns Hopkins University,} 2004.}
#'
#' @return A \eqn{d\times n} matrix that contains, column-wise, the sampled points from the convex polytope P.
#' @examples
#' # 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, n) {
    .Call(`_volesti_direct_sampling`, body, n)
}

#' Compute the exact volume of (a) a zonotope (b) an arbitrary simplex in V-representation or (c) if the volume is known and declared by the input object.
#'
#' Given a zonotope (as an object of class Zonotope), this function computes the sum of the absolute values of the determinants of all the \eqn{d \times d} submatrices of the \eqn{m\times d} matrix \eqn{G} that contains row-wise the \eqn{m} \eqn{d}-dimensional segments that define the zonotope.
#' For an arbitrary simplex that is given in V-representation this function computes the absolute value of the determinant formed by the simplex's points assuming it is shifted to the origin.
#'
#' @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(2, 5)
#' vol = exact_vol(Z)
#'
#' \donttest{# compute the exact volume of a 2-d arbitrary simplex
#' V = matrix(c(2,3,-1,7,0,0),ncol = 2, nrow = 3, byrow = TRUE)
#' P = Vpolytope(V = V)
#' vol = exact_vol(P)
#' }
#'
#' # compute the exact volume the 10-dimensional cross polytope
#' P = gen_cross(10,'V')
#' vol = exact_vol(P)
#' @export
exact_vol <- function(P) {
    .Call(`_volesti_exact_vol`, P)
}

#' 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 simplex.
#'
#' @param a A \eqn{d}-dimensional vector that defines the direction of the hyperplane.
#' @param z0 The scalar that defines the half-space.
#'
#' @references \cite{Varsi, Giulio,
#' \dQuote{The multidimensional content of the frustum of the simplex,} \emph{Pacific J. Math. 46, no. 1, 303--314,} 1973.}
#'
#' @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 simplex that is contained in the intersection of a given half-space and the simplex.
#'
#' @examples
#' # compute the frustum of H: -x1+x2<=0
#' a=c(-1,1)
#' z0=0
#' frustum = frustum_of_simplex(a, z0)
#' @export
frustum_of_simplex <- function(a, z0) {
    .Call(`_volesti_frustum_of_simplex`, 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{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 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.
#'
#' @examples
#' # compute the Chebychev ball of the 2d unit simplex
#' P = gen_simplex(2,'H')
#' ball_vec = inner_ball(P)
#'
#' # compute an inscribed ball of the 3-dimensional unit cube in V-representation
#' P = gen_cube(3, 'V')
#' ball_vec = inner_ball(P)
#' @export
inner_ball <- function(P) {
    .Call(`_volesti_inner_ball`, P)
}

#'  An internal Rccp function to read a SDPA format file
#'
#' @param input_file Name of the input file
#'
#' @keywords internal
#'
#' @return A list with two named items: an item "matrices" which is a list of the matrices and an vector "objFunction"
load_sdpa_format_file <- function(input_file = NULL) {
    .Call(`_volesti_load_sdpa_format_file`, input_file)
}

#' An internal Rccp function as a polytope generator
#'
#' @param kind_gen An integer to declare the type of the polytope.
#' @param Vpoly_gen A boolean parameter to declare if the requested polytope has to be in V-representation.
#' @param Zono_gen A boolean parameter to declare if the requested polytope has to be a zonotope.
#' @param dim_gen An integer to declare the dimension of the requested polytope.
#' @param m_gen An integer to declare the number of generators for the requested random zonotope or the number of vertices for a V-polytope.
#' @param seed Optional. A fixed seed for the random polytope generator.
#'
#' @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) {
    .Call(`_volesti_poly_gen`, kind_gen, Vpoly_gen, Zono_gen, dim_gen, m_gen, seed)
}

#'  An internal Rccp function for the random rotation of a convex polytope
#'
#' @param P A convex polytope (H-, V-polytope or a zonotope).
#' @param T Optional. A rotation matrix.
#' @param seed Optional. A fixed seed for the random linear map generator.
#'
#' @keywords internal
#'
#' @return A matrix that describes the rotated polytope
rotating <- function(P, T = NULL, seed = NULL) {
    .Call(`_volesti_rotating`, P, T, seed)
}

#' Internal rcpp function for the rounding of a convex polytope
#'
#' @param P A convex polytope (H- or V-representation or zonotope).
#' @param settings A list to set the random walk and its walk length
#' @param seed Optional. A fixed seed for the number generator.
#'
#' @keywords internal
#'
#' @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, settings = NULL, seed = NULL) {
    .Call(`_volesti_rounding`, P, settings, seed)
}

#' 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- 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) 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 per generated point for the random walk. The default value is 1.}
#' \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 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{\code{seed} }{ A fixed seed for the number generator.}
#' }
#' @param distribution Optional. A list that declares the target density and some related parameters as follows:
#' \itemize{
#' \item{\code{density} }{ A string: (a) \code{'uniform'} for the uniform distribution or b) \code{'gaussian'} for the multidimensional spherical distribution. The default target distribution is uniform.}
#' \item{\code{variance} }{ The variance of the multidimensional spherical gaussian. The default value is 1.}
#'  \item{\code{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()}.}
#' }
#'
#' @return A \eqn{d\times n} matrix that contains, column-wise, the sampled points from the convex polytope P.
#' @examples
#' # uniform distribution from the 3d unit cube in H-representation using ball walk
#' 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
#' A = matrix(c(-1,0,0,-1,1,1), ncol=2, nrow=3, byrow=TRUE)
#' b = c(0,0,1)
#' P = Hpolytope(A = A, b = b)
#' points = sample_points(P, n = 100, distribution = list("density" = "gaussian", "variance" = 2))
#'
#' # uniform points from the boundary of a 2-dimensional random H-polytope
#' P = gen_rand_hpoly(2,20)
#' points = sample_points(P, n = 100, random_walk = list("walk" = "BRDHR"))
#'
#' @export
sample_points <- function(P, n, random_walk = NULL, distribution = NULL) {
    .Call(`_volesti_sample_points`, P, n, random_walk, distribution)
}

#' 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 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 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{'CB'} for CB algorithm, b) \code{'SoB'} for SOB algorithm or b) \code{'CG'} for CG algorithm. The defalut 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 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{\code{seed} }{ A fixed seed for the number generator.}
#' }
#' @param rounding A boolean parameter for rounding. The default value is \code{FALSE}.
#'
#' @references \cite{I.Z.Emiris and V. Fisikopoulos,
#' \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.}
#'
#'
#' @return The approximation of the volume of a convex polytope.
#' @examples
#'
#' # 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" = 2))
#'
#' @export
volume <- function(P, settings = NULL, rounding = FALSE) {
    .Call(`_volesti_volume`, P, settings, rounding)
}

#' Write a SDPA format file
#'
#' Outputs a spectrahedron (the matrices defining a linear matrix inequality) and a vector (the objective function)
#' to a SDPA format file.
#'
#' @param spectrahedron A spectrahedron in n dimensions; must be an object of class Spectrahedron
#' @param objective_function A numerical vector of length n
#' @param output_file Name of the output file
#'
#' @examples
#' \donttest{
#' A0 = matrix(c(-1,0,0,0,-2,1,0,1,-2), nrow=3, ncol=3, byrow = TRUE)
#' A1 = matrix(c(-1,0,0,0,0,1,0,1,0), nrow=3, ncol=3, byrow = TRUE)
#' A2 = matrix(c(0,0,-1,0,0,0,-1,0,0), nrow=3, ncol=3, byrow = TRUE)
#' lmi = list(A0, A1, A2)
#' S = Spectrahedron(matrices = lmi)
#' objFunction = c(1,1)
#' write_sdpa_format_file(S, objFunction, "output.txt")
#' }
#' @export
write_sdpa_format_file <- function(spectrahedron, objective_function, output_file) {
    invisible(.Call(`_volesti_write_sdpa_format_file`, spectrahedron, objective_function, output_file))
}

#' An internal Rccp function for the over-approximation of a zonotope
#'
#' @param Z A zonotope.
#' @param fit_ratio Optional. A boolean parameter to request the computation of the ratio of fitness.
#' @param settings Optional. A list that declares the values of the parameters of CB algorithm.
#' @param seed Optional. A fixed seed for the number generator.
#'
#' @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) {
    .Call(`_volesti_zono_approx`, Z, fit_ratio, settings, seed)
}