diff --git a/DESCRIPTION b/DESCRIPTION index 646c892f..a91a2ee9 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -20,5 +20,5 @@ Imports: methods, stats LinkingTo: Rcpp, RcppEigen, BH Suggests: testthat Encoding: UTF-8 -RoxygenNote: 7.2.3 +RoxygenNote: 7.3.1 BugReports: https://github.com/GeomScale/volesti/issues diff --git a/NAMESPACE b/NAMESPACE index e12919ab..0d88effd 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -2,9 +2,13 @@ export(compute_indicators) export(copula) +export(dinvweibull_with_loc) export(direct_sampling) +export(ess) +export(estimtate_lipschitz_constant) export(exact_vol) export(frustum_of_simplex) +export(gen_birkhoff) export(gen_cross) export(gen_cube) export(gen_prod_simplex) @@ -13,12 +17,20 @@ export(gen_rand_vpoly) export(gen_rand_zonotope) export(gen_simplex) export(gen_skinny_cube) +export(geweke) export(inner_ball) +export(loadSdpaFormatFile) +export(ode_solve) +export(pinvweibull_with_loc) +export(psrf_multivariate) +export(psrf_univariate) +export(raftery) export(read_sdpa_format_file) export(rotate_polytope) export(round_polytope) export(sample_points) export(volume) +export(writeSdpaFormatFile) export(write_sdpa_format_file) export(zonotope_approximation) exportClasses(Hpolytope) diff --git a/R/RcppExports.R b/R/RcppExports.R index c3f192ae..a4680226 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -99,8 +99,8 @@ ess <- function(samples) { #' 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$new(V) +#' A = matrix(c(2,3,-1,7,0,0),ncol = 2, nrow = 3, byrow = TRUE) +#' P = Vpolytope(V = A) #' vol = exact_vol(P) #' } #' @@ -176,6 +176,17 @@ inner_ball <- function(P, lpsolve = NULL) { .Call(`_volesti_inner_ball`, P, lpsolve) } +#' 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) +} + #' Solve an ODE of the form dx^n / dt^n = F(x, t) #' #' @param n The number of steps. @@ -305,7 +316,7 @@ rounding <- function(P, method = NULL, seed = NULL) { #' \item{\code{BaW_rad} }{ The radius for the ball walk.} #' \item{\code{L} }{ The maximum length of the billiard trajectory or the radius for the step of dikin, vaidya or john walk.} #' \item{\code{solver} }{ Specify ODE solver for logconcave sampling. Options are i) leapfrog, ii) euler iii) runge-kutta iv) richardson} -#' \item{\code{step_size }{ Optionally chosen step size for logconcave sampling. Defaults to a theoretical value if not provided.} +#' \item{\code{step_size} }{ Optionally chosen step size for logconcave sampling. Defaults to a theoretical value if not provided.} #' } #' @param distribution Optional. A list that declares the target density and some related parameters as follows: #' \itemize{ @@ -347,7 +358,7 @@ rounding <- function(P, method = NULL, seed = NULL) { #' # gaussian distribution from the 2d unit simplex in H-representation with variance = 2 #' A = matrix(c(-1,0,0,-1,1,1), ncol=2, nrow=3, byrow=TRUE) #' b = c(0,0,1) -#' P = Hpolytope$new(A,b) +#' P = Hpolytope(A = A, b = b) #' points = sample_points(P, n = 100, distribution = list("density" = "gaussian", "variance" = 2)) #' #' # uniform points from the boundary of a 2-dimensional random H-polytope @@ -376,7 +387,7 @@ sample_points <- function(P, n, random_walk = NULL, distribution = NULL, seed = #' 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$new(lmi); +#' S = Spectrahedron(matrices = lmi) #' objFunction = c(1,1) #' writeSdpaFormatFile(S, objFunction, "output.txt") #' } @@ -412,6 +423,7 @@ loadSdpaFormatFile <- function(inputFile = NULL) { #' \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{250} for CB with BiW and \eqn{400+3d^2} for CB and any other random walk and \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 Optional. A string parameter to request a rounding method to be applied in the input polytope before volume computation: a) \code{'min_ellipsoid'}, b) \code{'svd'}, c) \code{'max_ellipsoid'} and d) \code{'none'} for no rounding. #' @param seed Optional. A fixed seed for the number generator. @@ -439,8 +451,32 @@ loadSdpaFormatFile <- function(inputFile = NULL) { #' pair_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) +volume <- function(P, settings = NULL, rounding = NULL) { + .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 diff --git a/R/compute_indicators.R b/R/compute_indicators.R index fe82cdee..695bb191 100644 --- a/R/compute_indicators.R +++ b/R/compute_indicators.R @@ -1,46 +1,75 @@ #' 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 \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. +#' 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. -#' @param m Optional. The number of slices for the copula. The default value is 100. -#' @param n Optional. The number of points to sample. The default value is \eqn{5\cdot 10^5}. -#' @param nwarning Optional. The number of consecutive indicators larger than 1 required to declare a warning period. The default value is 60. -#' @param ncrisis Optional. The number of consecutive indicators larger than 1 required to declare a crisis period. The default value is 100. -#' @param seed Optional. A fixed seed for the number generator. +#' @param parameters A list to set a parameterization. +#' \itemize{ +#' \item{win_length }{ The length of the sliding window. The default value is 60.} +#' \item{m } { The number of slices for the copula. The default value is 100.} +#' \item{n }{ The number of points to sample. The default value is \eqn{5\cdot 10^5}.} +#' \item{nwarning }{ The number of consecutive indicators larger than 1 required to declare a warning period. The default value is 60.} +#' \item{ncrisis }{ The number of consecutive indicators larger than 1 required to declare a crisis period. The default value is 100.} +#' \item{seed }{ 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 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 +#' @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) +#' market_states_and_indicators = compute_indicators(asset_returns, +#' parameters = list("win_length" = 10, "m" = 10, "n" = 10000, "nwarning" = 2, "ncrisis" = 3)) #' #' @export -compute_indicators <- function(returns, win_length = NULL, m = NULL, n = NULL, nwarning = NULL, ncrisis = NULL, seed = NULL) { - - if (is.null(win_length)) win_length = 60 - if (is.null(m)) m = 100 - if (is.null(n)) n = 500000 - if (is.null(nwarning)) nwarning = 60 - if (is.null(ncrisis)) ncrisis = 100 - +#' @useDynLib volesti, .registration=TRUE +#' @importFrom Rcpp evalCpp +#' @importFrom Rcpp loadModule +#' @importFrom "utils" "read.csv" +#' @importFrom "stats" "cov" +#' @importFrom "methods" "new" +compute_indicators <- function(returns, parameters = list("win_length" = 60, "m" = 100, "n" = 500000, "nwarning" = 60, "ncrisis" = 100)) { + + win_length = 60 + if (!is.null(parameters$win_length)) { + win_length = parameters$win_length + } + m=100 + if (!is.null(parameters$m)){ + m = parameters$m + } + n = 500000 + if (!is.null(parameters$n)){ + n = parameters$n + } + nwarning = 60 + if (!is.null(parameters$nwarning)) { + nwarning = parameters$nwarning + } + ncrisis = 100 + if (!is.null(parameters$ncrisis)) { + ncrisis = parameters$ncrisis + } + seed = NULL + if (!is.null(parameters$seed)) { + seed = parameters$seed + } + nrows = dim(returns)[1] nassets = dim(returns)[2] wl = win_length-1 - + indicators = c() for (i in 1:(nrows-wl)) { - + Win=i:(i+wl) E = cov(returns[Win,]) - + compRet = rep(1,nassets) for (j in 1:nassets) { for (k in Win) { @@ -48,11 +77,11 @@ compute_indicators <- function(returns, win_length = NULL, m = NULL, n = NULL, n } compRet[j] = compRet[j] - 1 } - + cop = copula(r1 = compRet, sigma = E, m = m, n = n, seed = seed) blue_mass = 0 red_mass = 0 - + for (row in 1:m) { for (col in 1:m) { if (row-col<=0.2*m && row-col>=-0.2*m) { @@ -76,7 +105,7 @@ compute_indicators <- function(returns, win_length = NULL, m = NULL, n = NULL, n col = rep("normal", N) for (i in 1:N) { - + if(indicators[i]>1 && !set_index){ index = i set_index = TRUE @@ -98,7 +127,7 @@ compute_indicators <- function(returns, win_length = NULL, m = NULL, n = NULL, n col[index:i] = "crisis" } } - + return(list("indicators" = indicators, market_states = col)) - -} + +} \ No newline at end of file diff --git a/R/file_to_polytope.R b/R/file_to_polytope.R deleted file mode 100644 index 19dd61ac..00000000 --- a/R/file_to_polytope.R +++ /dev/null @@ -1,113 +0,0 @@ -#' 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 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. -#' -#' @return 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. -#' -#' @export -#' @useDynLib volesti, .registration=TRUE -#' @importFrom Rcpp evalCpp -#' @importFrom Rcpp loadModule -#' @importFrom "utils" "read.csv" -#' @importFrom "stats" "cov" -#' @importFrom "methods" "new" -#' @exportPattern "^[[:alpha:]]+" -file_to_polytope <- function(path, zonotope = FALSE){ - - ineorext=substr(path, start = nchar(path) - 2, stop = nchar(path)) - if(ineorext!="ine" && ineorext!="ext") { - stop("Only ine or ext files can be handled by this function!") - } - P = read.csv(path) - r = as.character(P[3,1]) - count_sp = 1 - str = "" - beg = 0 - for (j in 1:nchar(r)) { - if (substr(r, start=j, stop=j) == " ") { - beg = beg + 1 - } else { - break - } - } - for (i in seq(from= beg + 1, to=nchar(r), by=1)) { - if (substr(r, start=i, stop=i) == " ") { - if (count_sp == 1) { - m = as.numeric(str) - str = "" - count_sp = count_sp + 1 - } else { - d = as.numeric(str) - str = "" - break - } - } else { - str = paste0(str, substr(r, start=i, stop=i)) - } - } - A = rep(0,d) - A[1] = m - A[2] = d - newrow = rep(0,d) - for (i in 4:(dim(P)[1] - 2)) { - r = P[i,1] - r = as.character(r) - str = "" - count = 1 - beg = 0 - for (j in 1:nchar(r)) { - if(substr(r, start=j, stop=j)==" "){ - beg = beg + 1 - } else { - break - } - } - sp_bef = FALSE - for (j in seq(from=beg + 1, to=nchar(r), by=1)) { - if (substr(r, start=j, stop=j) == " "){ - if (sp_bef) { - next - } - sp_bef = TRUE - newrow[count] = as.numeric(str) - str = "" - count = count + 1 - } else { - str = paste0(str, substr(r, start=j, stop=j)) - sp_bef = FALSE - if (j == nchar(r)) { - newrow[count] = as.numeric(str) - } - } - } - A = rbind(A,newrow) - newrow = rep(0,d) - } - A = matrix(A, ncol=dim(A)[2]) # now matrix A is in ine or ext format - - # remove first row - A = A[-c(1),] - - # first column is the vector b - b = A[,1] - - # remove first column - A2 = A[,-c(1)] - - if(ineorext=="ine") { - P = Hpolytope$new(-A2,b) - } else { - if(!missing(zonotope)){ - if(zonotope) { - P = Zonotope$new(A2) - return(P) - } - } - P = Vpolytope$new(A2) - } - return(P) -} diff --git a/R/gen_birkhoff.R b/R/gen_birkhoff.R index 935f637e..85333caf 100644 --- a/R/gen_birkhoff.R +++ b/R/gen_birkhoff.R @@ -1,28 +1,28 @@ #' Generator function for Birkhoff polytope -#' +#' #' This function can be used to generate the full dimensional \eqn{n}-Birkhoff polytope in H-representation. #' The dimension of the generated polytope is \eqn{(n-1)^2}. -#' +#' #' @param n The order of the Birkhoff polytope -#' +#' #' @return A polytope class representing the full dimensional \eqn{n}-Birkhoff polytope in H-representation. -#' @examples +#' @examples #' # generate the Birkhoff polytope of order 5 #' P = gen_birkhoff(5) #' @export gen_birkhoff <- function(n) { - + kind_gen = 7 m_gen = 0 - + Mat = poly_gen(kind_gen, FALSE, FALSE, n, m_gen) - + # first column is the vector b b = Mat[,1] Mat = Mat[,-c(1)] - - P = Hpolytope$new(Mat, b) - + + P = Hpolytope(A = Mat, b = b) + return(P) - + } diff --git a/R/gen_cross.R b/R/gen_cross.R index d782f399..e9fde1f8 100644 --- a/R/gen_cross.R +++ b/R/gen_cross.R @@ -1,41 +1,41 @@ #' Generator function for cross polytopes -#' +#' #' 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. -#' +#' @param representation A string to declare the representation. It has to be \code{'H'} for H-representation or \code{'V'} for V-representation. Default valus is 'H'. +#' #' @return A polytope class representing a cross polytope in H- or V-representation. -#' @examples +#' @examples #' # generate a 10-dimensional cross polytope in H-representation #' P = gen_cross(5, 'H') -#' +#' #' # generate a 15-dimension cross polytope in V-representation #' P = gen_cross(15, 'V') #' @export -gen_cross <- function(dimension, representation) { - +gen_cross <- function(dimension, representation = 'H') { + kind_gen = 2 m_gen = 0 - if (representation == "V") { + if (representation == 'V') { Vpoly_gen = TRUE - } else if (representation == "H") { + } else if (representation == 'H') { Vpoly_gen = FALSE } else { stop('Not a known representation.') } - + Mat = poly_gen(kind_gen, Vpoly_gen, FALSE, dimension, m_gen) - + # first column is the vector b - b = Mat[,1] - Mat = Mat[,-c(1)] - + b = Mat[, 1] + Mat = Mat[, -c(1), drop = FALSE] + if (Vpoly_gen) { - P = Vpolytope$new(Mat, 2^dimension / prod(1:dimension)) + P = Vpolytope(V = Mat, volume = 2^dimension / prod(1:dimension)) } else { - P = Hpolytope$new(-Mat, b, 2^dimension / prod(1:dimension)) + P = Hpolytope(A = -Mat, b = b, volume = 2^dimension / prod(1:dimension)) } - + return(P) -} +} \ No newline at end of file diff --git a/R/gen_cube.R b/R/gen_cube.R index bcbf1f02..f189b720 100644 --- a/R/gen_cube.R +++ b/R/gen_cube.R @@ -1,41 +1,40 @@ #' Generator function for hypercubes -#' +#' #' 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. -#' +#' @param representation A string to declare the representation. It has to be \code{'H'} for H-representation or \code{'V'} for V-representation. Default valus is 'H'. +#' #' @return A polytope class representing the unit \eqn{d}-dimensional hypercube in H- or V-representation. -#' @examples +#' @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(5, 'V') #' @export -gen_cube <- function(dimension, representation) { - +gen_cube <- function(dimension, representation = 'H') { + kind_gen = 1 m_gen = 0 - if (representation == "V") { + if (representation == 'V') { Vpoly_gen = TRUE - } else if (representation == "H") { + } else if (representation == 'H') { Vpoly_gen = FALSE } else { stop('Not a known representation.') } - + Mat = poly_gen(kind_gen, Vpoly_gen, FALSE, dimension, m_gen) - + # first column is the vector b - b = Mat[,1] - Mat = Mat[,-c(1)] + b = Mat[, 1] + Mat = Mat[, -c(1), drop = FALSE] if (Vpoly_gen) { - P = Vpolytope$new(Mat, 2^dimension) + P = Vpolytope(V = Mat, volume = 2^dimension) } else { - P = Hpolytope$new(-Mat, b, 2^dimension) + P = Hpolytope(A = -Mat, b = b, volume = 2^dimension) } - + return(P) - } diff --git a/R/gen_prod_simplex.R b/R/gen_prod_simplex.R index a0505880..3f724e8a 100644 --- a/R/gen_prod_simplex.R +++ b/R/gen_prod_simplex.R @@ -1,29 +1,29 @@ #' Generator function for product of simplices -#' +#' #' 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. -#' +#' #' @return A polytope class representing the product of the two \eqn{d}-dimensional unit simplices in H-representation. -#' +#' #' @examples #' # generate a product of two 5-dimensional simplices. #' P = gen_prod_simplex(5) #' @export gen_prod_simplex <- function(dimension) { - + kind_gen = 4 m_gen = 0 Vpoly_gen = FALSE - + Mat = poly_gen(kind_gen, Vpoly_gen, FALSE, dimension, m_gen) # first column is the vector b - b = Mat[,1] - Mat = Mat[,-c(1)] - - P = Hpolytope$new(-Mat, b, (1/prod(1:dimension))^2) - + b = Mat[, 1] + Mat = Mat[, -c(1), drop = FALSE] + + P = Hpolytope(A = -Mat, b = b, volume = (1/prod(1:dimension))^2) + return(P) - -} + +} \ No newline at end of file diff --git a/R/gen_rand_hpoly.R b/R/gen_rand_hpoly.R index 2785a28b..a7208b82 100644 --- a/R/gen_rand_hpoly.R +++ b/R/gen_rand_hpoly.R @@ -1,28 +1,45 @@ #' Generator function for random H-polytopes -#' +#' #' 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. -#' @param seed Optional. A fixed seed for the generator. -#' +#' @param generator A list that could contain two elements. +#' \itemize{ +#' \item{constants }{ To declare how to set the constants \eqn{b_i} for each facets: (i) 'sphere', each hyperplane is tangent to the hypersphere of radius 10, (ii) 'uniform' for each \eqn{b_i} the generator picks a uniform number from \eqn{(0,1)}. The defalut value is 'sphere'.} +#' \item{seed }{ Optional. A fixed seed for the number generator.} +#' } +#' #' @return A polytope class representing a H-polytope. -#' @examples +#' @examples #' # generate a 10-dimensional polytope with 50 facets #' P = gen_rand_hpoly(10, 50) #' @export -gen_rand_hpoly <- function(dimension, nfacets, seed = NULL) { - - kind_gen = 6 +gen_rand_hpoly <- function(dimension, nfacets, generator = list('constants' = 'sphere')) { + + seed = NULL + if (!is.null(generator$seed)) { + seed = generator$seed + } + + if (is.null(generator$constants)) { + kind_gen = 6 + } else if (generator$constants == 'sphere'){ + kind_gen = 6 + } else if (generator$constants == 'uniform') { + kind_gen = 7 + } else { + stop("Wrong generator!") + } Vpoly_gen = FALSE - + Mat = poly_gen(kind_gen, Vpoly_gen, FALSE, dimension, nfacets, seed) - + # first column is the vector b - b = Mat[,1] - Mat = Mat[,-c(1)] - - P = Hpolytope$new(Mat, b) - + b = Mat[, 1] + Mat = Mat[, -c(1), drop = FALSE] + + P = Hpolytope(A = Mat, b = b) + return(P) -} +} \ No newline at end of file diff --git a/R/gen_rand_vpoly.R b/R/gen_rand_vpoly.R index 31cc561c..e67f93aa 100644 --- a/R/gen_rand_vpoly.R +++ b/R/gen_rand_vpoly.R @@ -1,36 +1,44 @@ #' Generator function for random V-polytopes -#' +#' #' 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. -#' @param generator The body that the generator samples uniformly the vertices from: (a) 'cube' or (b) 'sphere'. -#' @param seed Optional. A fixed seed for the generator. -#' +#' @param generator A list that could contain two elements. +#' \itemize{ +#' \item{body }{ the body that the generator samples uniformly the vertices from: (i) 'cube' or (ii) 'sphere', the default value is 'sphere'.} +#' \item{seed }{ Optional. A fixed seed for the number generator.} +#' } +#' #' @return A polytope class representing a V-polytope. -#' @examples +#' @examples #' # generate a 10-dimensional polytope defined as the convex hull of 25 random vertices #' P = gen_rand_vpoly(10, 25) #' @export -gen_rand_vpoly <- function(dimension, nvertices, generator = NULL, seed = NULL) { - - kind_gen = 4 - - if(!is.null(generator)){ - if (generator == 'cube'){ - kind_gen = 5 - } else if (generator != 'sphere') { - stop("Wrong generator!") - } +gen_rand_vpoly <- function(dimension, nvertices, generator = list('body' = 'sphere')) { + + seed = NULL + if (!is.null(generator$seed)) { + seed = generator$seed + } + + if (is.null(generator$body)) { + kind_gen = 4 + } else if (generator$body == 'cube'){ + kind_gen = 5 + } else if (generator$body == 'sphere') { + kind_gen = 4 + } else { + stop("Wrong generator!") } - + Mat = poly_gen(kind_gen, TRUE, FALSE, dimension, nvertices, seed) # first column is the vector b - b = Mat[,1] - Mat = Mat[,-c(1)] - - P = Vpolytope$new(Mat) - + b = Mat[, 1] + Mat = Mat[, -c(1), drop = FALSE] + + P = Vpolytope(V = Mat) + return(P) -} +} \ No newline at end of file diff --git a/R/gen_rand_zonotope.R b/R/gen_rand_zonotope.R index 3a2c8c1b..a7286dc8 100644 --- a/R/gen_rand_zonotope.R +++ b/R/gen_rand_zonotope.R @@ -1,40 +1,48 @@ #' Generator function for zonotopes -#' +#' #' 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. -#' @param generator The distribution to pick the length of each segment from \eqn{[0,100]}: (a) 'uniform', (b) 'gaussian' or (c) 'exponential'. -#' @param seed Optional. A fixed seed for the generator. -#' +#' @param generator A list that could contain two elements. +#' \itemize{ +#' \item{distribution }{ the distribution to pick the length of each segment from \eqn{[0,100]}: (i) 'uniform', (ii) 'gaussian' or (iii) 'exponential', the default value is 'uniform.} +#' \item {seed }{ Optional. A fixed seed for the number generator.} +#' } +#' #' @return A polytope class representing a zonotope. #' -#' @examples +#' @examples #' # generate a 10-dimensional zonotope defined by the Minkowski sum of 20 segments #' P = gen_rand_zonotope(10, 20) #' @export -gen_rand_zonotope <- function(dimension, nsegments, generator = NULL, seed = NULL) { - - kind_gen = 1 - - if (!is.null(generator)) { - if (generator == 'gaussian') { - kind_gen = 2 - } else if (generator == 'exponential') { - kind_gen = 3 - } else if (generator != 'uniform'){ - stop("Wrong generator!") - } +gen_rand_zonotope <- function(dimension, nsegments, generator = list('distribution' = 'uniform')) { + + seed = NULL + if (!is.null(generator$seed)) { + seed = generator$seed } - + + if (is.null(generator$distribution)) { + kind_gen = 1 + } else if (generator$distribution == 'gaussian') { + kind_gen = 2 + } else if (generator$distribution == 'exponential') { + kind_gen = 3 + } else if (generator$distribution == 'uniform'){ + kind_gen = 1 + } else { + stop("Wrong generator!") + } + Mat = poly_gen(kind_gen, FALSE, TRUE, dimension, nsegments, seed) - + # first column is the vector b - b = Mat[,1] - Mat = Mat[,-c(1)] - - P = Zonotope$new(Mat) + b = Mat[, 1] + Mat = Mat[, -c(1), drop = FALSE] + + P = Zonotope(G = Mat) return(P) -} +} \ No newline at end of file diff --git a/R/gen_simplex.R b/R/gen_simplex.R index aa2ca4a6..99fe29d5 100644 --- a/R/gen_simplex.R +++ b/R/gen_simplex.R @@ -1,20 +1,20 @@ #' Generator function for simplices -#' +#' #' 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. -#' +#' @param representation A string to declare the representation. It has to be \code{'H'} for H-representation or \code{'V'} for V-representation. Default valus is 'H'. +#' #' @return A polytope class representing the \eqn{d}-dimensional unit simplex in H- or V-representation. #' @examples #' # generate a 10-dimensional simplex in H-representation #' PolyList = gen_simplex(10, 'H') -#' +#' #' # generate a 20-dimensional simplex in V-representation #' P = gen_simplex(20, 'V') #' @export -gen_simplex <- function(dimension, representation) { - +gen_simplex <- function(dimension, representation = 'H') { + kind_gen = 3 m_gen = 0 if (representation == "V") { @@ -24,18 +24,18 @@ gen_simplex <- function(dimension, representation) { } else { stop('Not a known representation.') } - + Mat = poly_gen(kind_gen, Vpoly_gen, FALSE, dimension, m_gen) # first column is the vector b - b = Mat[,1] - Mat = Mat[,-c(1)] - + b = Mat[, 1] + Mat = Mat[, -c(1), drop = FALSE] + if (Vpoly_gen) { - P = Vpolytope$new(Mat, 1/prod(1:dimension)) + P = Vpolytope(V = Mat, volume = 1/prod(1:dimension)) } else { - P = Hpolytope$new(-Mat, b, 1/prod(1:dimension)) + P = Hpolytope(A = -Mat, b = b, volume = 1/prod(1:dimension)) } - + return(P) -} +} \ No newline at end of file diff --git a/R/gen_skinny_cube.R b/R/gen_skinny_cube.R index 0a1df0eb..2dde1d1e 100644 --- a/R/gen_skinny_cube.R +++ b/R/gen_skinny_cube.R @@ -1,9 +1,9 @@ #' Generator function for skinny hypercubes -#' +#' #' 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. -#' +#' #' @return A polytope class representing the \eqn{d}-dimensional skinny hypercube in H-representation. #' #' @examples @@ -11,18 +11,18 @@ #' P = gen_skinny_cube(10) #' @export gen_skinny_cube <- function(dimension) { - + kind_gen = 5 m_gen = 0 Vpoly_gen = FALSE - + Mat = poly_gen(kind_gen, Vpoly_gen, FALSE, dimension, m_gen) # first column is the vector b - b = Mat[,1] - Mat = Mat[,-c(1)] - - P = Hpolytope$new(-Mat, b, 2^(dimension -1)*200) - + b = Mat[, 1] + Mat = Mat[, -c(1), drop = FALSE] + + P = Hpolytope(A = -Mat, b = b, volume = 2^(dimension -1)*200) + return(P) -} +} \ No newline at end of file diff --git a/R/read_sdpa_file.R b/R/read_sdpa_file.R index 83c10c52..4c78523a 100644 --- a/R/read_sdpa_file.R +++ b/R/read_sdpa_file.R @@ -9,7 +9,7 @@ #' #' @examples #' path = system.file('extdata', package = 'volesti') -#' l = readSdpaFormatFile(paste0(path,'/sdpa_n2m3.txt')) +#' l = read_sdpa_format_file(paste0(path,'/sdpa_n2m3.txt')) #' Spectrahedron = l$spectrahedron #' objFunction = l$objFunction #' @export @@ -18,9 +18,9 @@ #' @importFrom Rcpp loadModule #' @importFrom "methods" "new" #' @exportPattern "^[[:alpha:]]+" -readSdpaFormatFile <- function(path){ - l = loadSdpaFormatFile(path) - S = Spectrahedron$new(l$matrices) +read_sdpa_format_file <- function(path){ + l = load_sdpa_format_file(path) + S = Spectrahedron(matrices = l$matrices) return(list("spectrahedron"=S, "objFunction"= l$objFunction)) -} +} \ No newline at end of file diff --git a/R/rotate_polytope.R b/R/rotate_polytope.R index 12ce246b..b901ee14 100644 --- a/R/rotate_polytope.R +++ b/R/rotate_polytope.R @@ -1,14 +1,13 @@ #' Apply a random rotation to a convex polytope (H-polytope, V-polytope, zonotope or intersection of two V-polytopes) -#' +#' #' Given a convex H- or V- polytope or a zonotope or an intersection of two V-polytopes as input, this function applies (a) a random rotation or (b) a given rotation by an input matrix \eqn{T}. -#' +#' #' @param P A convex polytope. It is an object from class (a) Hpolytope, (b) Vpolytope, (c) Zonotope, (d) intersection of two V-polytopes. -#' @param T Optional. A \eqn{d\times d} rotation matrix. -#' @param seed Optional. A fixed seed for the random linear map generator. -#' +#' @param rotation A list that contains (a) the rotation matrix T and (b) the 'seed' to set a spesific seed for the number generator. +#' #' @return A list that contains the rotated polytope and the matrix \eqn{T} of the linear transformation. #' -#' @details Let \eqn{P} be the given polytope and \eqn{Q} the rotated one and \eqn{T} be the matrix of the linear transformation. +#' @details Let \eqn{P} be the given polytope and \eqn{Q} the rotated one and \eqn{T} be the matrix of the linear transformation. #' \itemize{ #' \item{If \eqn{P} is in H-representation and \eqn{A} is the matrix that contains the normal vectors of the facets of \eqn{Q} then \eqn{AT} contains the normal vactors of the facets of \eqn{P}.} #' \item{If \eqn{P} is in V-representation and \eqn{V} is the matrix that contains column-wise the vertices of \eqn{Q} then \eqn{T^TV} contains the vertices of \eqn{P}.} @@ -17,41 +16,60 @@ #' } #' @examples #' # rotate a H-polytope (2d unit simplex) -#' P = gen_simplex(2,'H') +#' P = gen_simplex(2, 'H') #' poly_matrix_list = rotate_polytope(P) -#' +#' #' # rotate a V-polytope (3d cube) #' P = gen_cube(3, 'V') #' poly_matrix_list = rotate_polytope(P) -#' +#' #' # rotate a 5-dimensional zonotope defined by the Minkowski sum of 15 segments -#' Z = gen_rand_zonotope(3,6) +#' Z = gen_rand_zonotope(3, 6) #' poly_matrix_list = rotate_polytope(Z) #' @export -rotate_polytope <- function(P, T = NULL, seed = NULL){ - +rotate_polytope <- function(P, rotation = list()) { + + seed = NULL + if (!is.null(rotation$seed)) { + seed = rotation$seed + } + + if (is.null(rotation$T)) { + T = NULL + } else { + T = rotation$T + } + #call rcpp rotating function Mat = rotating(P, T, seed) - - n = P$dimension - m=dim(Mat)[2]-n - Tr = Mat[,-c(1:(dim(Mat)[2]-n))] + + type = P@type + + if (type == 'Vpolytope') { + n = dim(P@V)[2] + }else if (type == 'Zonotope') { + n = dim(P@G)[2] + } else { + n = dim(P@A)[2] + } + + m = dim(Mat)[2] - n + Tr = Mat[, -c(1:(dim(Mat)[2]-n)), drop = FALSE] Tr = Tr[1:n, 1:n] - Mat = t(Mat[,1:m]) - + Mat = t(Mat[, 1:m]) + # first column is the vector b - b = Mat[,1] - + b = Mat[, 1] + # remove first column - A = Mat[,-c(1)] - - type = P$type - if (type == 2) { - PP = Vpolytope$new(A) - }else if (type == 3) { - PP = Zonotope$new(A) + A = Mat[, -c(1), drop = FALSE] + + if (type == 'Vpolytope') { + PP = Vpolytope(V = A) + }else if (type == 'Zonotope') { + PP = Zonotope(G = A) } else { - PP = Hpolytope$new(A, b) + PP = Hpolytope(A = A, b = b) } return(list("P" = PP, "T" = Tr)) -} +} \ No newline at end of file diff --git a/R/round_polytope.R b/R/round_polytope.R index f95649c6..e50ab7b9 100644 --- a/R/round_polytope.R +++ b/R/round_polytope.R @@ -1,11 +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 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 method Optional. The method to use for rounding, a) \code{'min_ellipsoid'} for the method based on mimimmum volume enclosing ellipsoid of a dataset, b) \code{'max_ellipsoid'} for the method based on maximum volume enclosed ellipsoid, (c) \code{'isotropy'} for the method based on svd decomposition. The default method is \code{'mee'} for all the representations. -#' @param seed Optional. A fixed seed for the number generator. -#' +#' @param settings Optional. A list of settings. +#' \itemize{ +#' \item{\code{method} }{ The method to use for rounding, a) \code{'min_ellipsoid'} for the method based on mimimmum volume enclosing ellipsoid of a dataset, b) \code{'max_ellipsoid'} for the method based on maximum volume enclosed ellipsoid, (c) \code{'isotropy'} for the method based on svd decomposition. The default method is \code{'mee'} for all the representations.} +#' \item{\code{seed} }{ Optional. A fixed seed for the number generator.} +#' } +#' #' @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, @@ -16,45 +19,48 @@ #' \dQuote{A practical volume algorithm,} \emph{Math. Prog. Comp.,} 2016.}, #' @references \cite{Yin Zhang and Liyan Gao, #' \dQuote{On Numerical Solution of the Maximum Volume Ellipsoid Problem,} \emph{SIAM Journal on Optimization,} 2003.}, -#' +#' #' #' @examples #' # round a 5d skinny cube #' P = gen_skinny_cube(5) #' listHpoly = round_polytope(P) -#' +#' #' # round a V-polytope (3d unit cube) #' P = gen_cube(3, 'V') #' ListVpoly = round_polytope(P) -#' +#' #' # round a 2-dimensional zonotope defined by 6 generators #' Z = gen_rand_zonotope(2,6) #' ListZono = round_polytope(Z) #' @export -round_polytope <- function(P, method = NULL, seed = NULL){ - - ret_list = rounding(P, method, seed) - +round_polytope <- function(P, settings = list()){ + + seed = NULL + if (!is.null(settings$seed)) { + seed = settings$seed + } + + ret_list = rounding(P, settings$method, seed) + #get the matrix that describes the polytope Mat = ret_list$Mat - + # first column is the vector b - b = Mat[,1] - + b = Mat[, 1] + # remove first column - A = Mat[,-c(1)] - - type = P$type - if (type == 2) { - PP = list("P" = Vpolytope$new(A), "T" = ret_list$T, "shift" = ret_list$shift, "round_value" = ret_list$round_value) - }else if (type == 3) { - PP = list("P" = Zonotope$new(A), "T" = ret_list$T, "shift" = ret_list$shift, "round_value" = ret_list$round_value) + A = Mat[, -c(1), drop = FALSE] + + type = P@type + if (type == 'Vpolytope') { + PP = list("P" = Vpolytope(V = A), "T" = ret_list$T, "shift" = ret_list$shift, "round_value" = ret_list$round_value) + }else if (type == 'Zonotope') { + PP = list("P" = Zonotope(G = A), "T" = ret_list$T, "shift" = ret_list$shift, "round_value" = ret_list$round_value) } else { - if (dim(P$Aeq)[1] > 0){ - PP = list("P" = Hpolytope$new(A,b), "T" = ret_list$T, "shift" = ret_list$shift, "round_value" = ret_list$round_value, "N" = ret_list$N, "N_shift" = ret_list$N_shift, "svd_prod" = ret_list$svd_prod) - } else { - PP = list("P" = Hpolytope$new(A,b), "T" = ret_list$T, "shift" = ret_list$shift, "round_value" = ret_list$round_value) - } + #if (dim(P$Aeq)[1] > 0){ + # PP = list("P" = Hpolytope(A = A, b = b), "T" = ret_list$T, "shift" = ret_list$shift, "round_value" = ret_list$round_value, "N" = ret_list$N, "N_shift" = ret_list$N_shift, "svd_prod" = ret_list$svd_prod) + PP = list("P" = Hpolytope(A = A, b = b), "T" = ret_list$T, "shift" = ret_list$shift, "round_value" = ret_list$round_value) } return(PP) } diff --git a/R/zonotope_approximation.R b/R/zonotope_approximation.R index 8c56f6d2..efbe73dc 100644 --- a/R/zonotope_approximation.R +++ b/R/zonotope_approximation.R @@ -1,42 +1,47 @@ #' 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. 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{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{win_len} }{ The length of the sliding window for CB algorithm. The default value is \eqn{250}.} #' \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{\code{seed} }{ Optional. A fixed seed for the number generator.} #' } -#' @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 -#' +#' #' @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,12) +#' Z = gen_rand_zonotope(2, 10) #' retList = zonotope_approximation(Z = Z) -#' +#' #' @export -zonotope_approximation <- function(Z, fit_ratio = NULL, settings = NULL, seed = NULL){ - +zonotope_approximation <- function(Z, fit_ratio = FALSE, settings = list('error' = 0.1, 'walk_length' = 1, 'win_len' = 250, 'hpoly' = FALSE)){ + + seed = NULL + if (!is.null(settings$seed)) { + seed = settings$seed + } + ret_list = zono_approx(Z, fit_ratio, settings, seed) - + Mat = ret_list$Mat - + # first column is the vector b - b = Mat[,1] - + b = Mat[, 1] + # remove first column - A = Mat[,-c(1)] - PP = list("P" = Hpolytope$new(A,b), "fit_ratio" = ret_list$fit_ratio) - + A = Mat[, -c(1), drop = FALSE] + PP = list("P" = Hpolytope(A = A, b = b), "fit_ratio" = ret_list$fit_ratio) + return(PP) - + } diff --git a/R/zzz.R b/R/zzz.R index 564a20b9..4bea4292 100644 --- a/R/zzz.R +++ b/R/zzz.R @@ -7,5 +7,5 @@ ## For R 2.15.1 and later this also works. Note that calling loadModule() triggers ## a load action, so this does not have to be placed in .onLoad() or evalqOnLoad(). -loadModule("polytopes", TRUE) -loadModule("spectrahedron", TRUE) +#loadModule("polytopes", TRUE) +#loadModule("spectrahedron", TRUE) diff --git a/man/Hpolytope-class.Rd b/man/Hpolytope-class.Rd new file mode 100644 index 00000000..1c5b2472 --- /dev/null +++ b/man/Hpolytope-class.Rd @@ -0,0 +1,27 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/HpolytopeClass.R +\docType{class} +\name{Hpolytope-class} +\alias{Hpolytope-class} +\alias{Hpolytope} +\title{An R class to represent an H-polytope} +\description{ +An H-polytope is a convex polytope defined by a set of linear inequalities or equivalently a \eqn{d}-dimensional H-polytope with \eqn{m} facets is defined by a \eqn{m\times d} matrix A and a \eqn{m}-dimensional vector b, s.t.: \eqn{Ax\leq b}. +} +\details{ +\describe{ + \item{A}{An \eqn{m\times d} numerical matrix.} + + \item{b}{An \eqn{m}-dimensional vector b.} + + \item{volume}{The volume of the polytope if it is known, \eqn{NaN} otherwise by default.} + + \item{type}{A character with default value 'Hpolytope', to declare the representation of the polytope.} +} +} +\examples{ +A = matrix(c(-1,0,0,-1,1,1), ncol=2, nrow=3, byrow=TRUE) +b = c(0,0,1) +P = Hpolytope(A = A, b = b) + +} diff --git a/man/Hpolytope.Rd b/man/Hpolytope.Rd deleted file mode 100644 index c0b82674..00000000 --- a/man/Hpolytope.Rd +++ /dev/null @@ -1,19 +0,0 @@ -\name{Hpolytope} -\alias{Hpolytope} -\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{P=\{x\ |\ Ax\leq b\} }. -} -\section{Fields}{ -\itemize{ -\item{\code{A} }{ A \eqn{m\times d} numerical matrix A} - -\item{\code{b} }{ \eqn{m}-dimensional vector b} - -\item{\code{type} }{ An integer that declares the representation of the polytope. For H-representation the default value is 1.} - -\item{\code{dimension} }{ The dimension of the polytope.} - -\item{\code{volume} }{ The volume of the polytope, if it is known.} -}} diff --git a/man/Rcpp_Hpolytope.Rd b/man/Rcpp_Hpolytope.Rd deleted file mode 100644 index 88fa5b76..00000000 --- a/man/Rcpp_Hpolytope.Rd +++ /dev/null @@ -1,28 +0,0 @@ -\docType{class} -\name{Rcpp_Hpolytope} -\alias{Rcpp_Hpolytope-class} -\alias{[,Rcpp_Hpolytope-method} -\alias{[,Rcpp_Hpolytope,ANY,ANY,ANY-method} -\alias{$<-,Rcpp_Hpolytope-method} -\alias{$,Rcpp_Hpolytope-method} -\alias{filepaths<-,Rcpp_Hpolytope-method} -\title{ -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{P=\{ Ax\leq b \} }. -} -\details{ -\describe{ -\item{\code{A} }{ A \eqn{m\times d} numerical matrix A} - -\item{\code{b} }{ \eqn{m}-dimensional vector b} - -\item{\code{type} }{ An integer that declares the representation of the polytope. For H-representation the default value is 1.} - -\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/man/Rcpp_Spectrahedron.Rd b/man/Rcpp_Spectrahedron.Rd deleted file mode 100644 index 3ac432e5..00000000 --- a/man/Rcpp_Spectrahedron.Rd +++ /dev/null @@ -1,21 +0,0 @@ -\docType{class} -\name{Rcpp_Spectrahedron} -\alias{Rcpp_Spectrahedron-class} -\alias{[,Rcpp_Spectrahedron-method} -\alias{[,Rcpp_Spectrahedron,ANY,ANY,ANY-method} -\alias{$<-,Rcpp_Spectrahedron-method} -\alias{$,Rcpp_Spectrahedron-method} -\alias{filepaths<-,Rcpp_Spectrahedron-method} -\title{ -An \code{Rcpp} class to represent spectrahedra, exposed to \code{R} via modules. -} -\description{ -A spectrahedron is a convex body defined by a linear matrix inequality of the form \eqn{A_0 + x_1 A_1 + ... + x_n A_n \preceq 0}. -The matrices \eqn{A_i} are symmetric \eqn{m \times m} real matrices and \eqn{\preceq 0} denoted negative semidefiniteness. -} -\details{ -\describe{ -\item{\code{matrices} }{A list with the matrices \eqn{A_0, A_1, ..., A_n}} - } -} - diff --git a/man/Rcpp_Vpolytope.Rd b/man/Rcpp_Vpolytope.Rd deleted file mode 100644 index 95c46c09..00000000 --- a/man/Rcpp_Vpolytope.Rd +++ /dev/null @@ -1,26 +0,0 @@ -\docType{class} -\name{Rcpp_Vpolytope} -\alias{Rcpp_Vpolytope-class} -\alias{[,Rcpp_Vpolytope-method} -\alias{[,Rcpp_Vpolytope,ANY,ANY,ANY-method} -\alias{$<-,Rcpp_Vpolytope-method} -\alias{$,Rcpp_Vpolytope-method} -\alias{filepaths<-,Rcpp_Vpolytope-method} -\title{ -An \code{Rcpp} class to represent V-polytopes, exposed to \code{R} via modules. -} -\description{ -A V-polytope is defined as the convex hull of \eqn{m} \eqn{d}-dimensional points which corresponds to its vertices. -} -\details{ -\describe{ - \item{\code{V} }{ A \eqn{m\times d} numerical matrix that contains the vertices row-wise} - - \item{\code{type} }{ An integer that declares the representation of the polytope. For V-representation the default value is 2.} - - \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/man/Rcpp_VpolytopeIntersection.Rd b/man/Rcpp_VpolytopeIntersection.Rd deleted file mode 100644 index ec308c3b..00000000 --- a/man/Rcpp_VpolytopeIntersection.Rd +++ /dev/null @@ -1,28 +0,0 @@ -\docType{class} -\name{Rcpp_VpolytopeIntersection} -\alias{Rcpp_VpolytopeIntersection-class} -\alias{[,Rcpp_VpolytopeIntersection-method} -\alias{[,Rcpp_VpolytopeIntersection,ANY,ANY,ANY-method} -\alias{$<-,Rcpp_VpolytopeIntersection-method} -\alias{$,Rcpp_VpolytopeIntersection-method} -\alias{filepaths<-,Rcpp_VpolytopeIntersection-method} -\title{ -An \code{Rcpp} class to represent the intersection of two V-polytope, exposed to \code{R} via modules. -} -\description{ -An intersection of two V-polytopes, \eqn{P_1}, \eqn{P_2}, is defined by the intersection of the two coresponding convex hulls. -} -\details{ -\describe{ -\item{\code{V1} }{ The numerical matrix that contains the vertices of \eqn{P_1} row-wise.} - -\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{dimension} }{ The dimension of the polytope.} - -\item{\code{volume} }{ The volume of the polytope, if it is known.} - } -} -\keyword{internal} diff --git a/man/Rcpp_Zonotope.Rd b/man/Rcpp_Zonotope.Rd deleted file mode 100644 index 8c010851..00000000 --- a/man/Rcpp_Zonotope.Rd +++ /dev/null @@ -1,26 +0,0 @@ -\docType{class} -\name{Rcpp_Zonotope} -\alias{Rcpp_Zonotope-class} -\alias{[,Rcpp_Zonotope-method} -\alias{[,Rcpp_Zonotope,ANY,ANY,ANY-method} -\alias{$<-,Rcpp_Zonotope-method} -\alias{$,Rcpp_Zonotope-method} -\alias{filepaths<-,Rcpp_Zonotope-method} -\title{ -An \code{Rcpp} class to represent zonotopes, exposed to \code{R} via modules. -} -\description{ -A zonotope is a convex polytope defined by the Minkowski sum of \eqn{m} \eqn{d}-dimensional segments. -} -\details{ -\describe{ - \item{\code{G} }{ A \eqn{m\times d} numerical matrix that contains the segments (or generators) row-wise} - -\item{\code{type} }{ An integer that declares the representation of the polytope. For zonotopes the default value is 3.} - -\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/man/Spectrahedron-class.Rd b/man/Spectrahedron-class.Rd new file mode 100644 index 00000000..2cc381a2 --- /dev/null +++ b/man/Spectrahedron-class.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/SpectrahedronClass.R +\docType{class} +\name{Spectrahedron-class} +\alias{Spectrahedron-class} +\alias{Spectrahedron} +\title{An R class to represent a Spectrahedron} +\description{ +A spectrahedron is a convex body defined by a linear matrix inequality of the form \eqn{A_0 + x_1 A_1 + ... + x_n A_n \preceq 0}. +The matrices \eqn{A_i} are symmetric \eqn{m \times m} real matrices and \eqn{\preceq 0} denoted negative semidefiniteness. +} +\details{ +\describe{ + \item{matrices}{A List that contains the matrices \eqn{A_0, A_1, ..., A_n}.} +} +} +\examples{ +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); + +} diff --git a/man/Spectrahedron.Rd b/man/Spectrahedron.Rd deleted file mode 100644 index 6c6df1d8..00000000 --- a/man/Spectrahedron.Rd +++ /dev/null @@ -1,12 +0,0 @@ -\name{Spectrahedron} -\alias{Spectrahedron} -\title{An \code{R} class to represent spectrahedra.} - -\description{ -A spectrahedron is a convex body defined by a linear matrix inequality of the form \eqn{A_0 + x_1 A_1 + ... + x_n A_n \preceq 0}. -The matrices \eqn{A_i} are symmetric \eqn{m \times m} real matrices and \eqn{\preceq 0} denoted negative semidefiniteness. -} -\section{Fields}{ -\itemize{ -\item{\code{matrices} }{A list with the matrices \eqn{A_0, A_1, ..., A_n}} -}} diff --git a/man/Vpolytope-class.Rd b/man/Vpolytope-class.Rd new file mode 100644 index 00000000..83fcc197 --- /dev/null +++ b/man/Vpolytope-class.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/VpolytopeClass.R +\docType{class} +\name{Vpolytope-class} +\alias{Vpolytope-class} +\alias{Vpolytope} +\title{An R class to represent a V-polytope} +\description{ +A V-polytope is a convex polytope defined by the set of its vertices. +} +\details{ +\describe{ + \item{V}{An \eqn{m\times d} numerical matrix that contains the vertices row-wise.} + + \item{volume}{The volume of the polytope if it is known, \eqn{NaN} otherwise by default.} + + \item{type}{A character with default value 'Vpolytope', to declare the representation of the polytope.} +} +} +\examples{ +V = matrix(c(2,3,-1,7,0,0),ncol = 2, nrow = 3, byrow = TRUE) +P = Vpolytope(V = V) + +} diff --git a/man/Vpolytope.Rd b/man/Vpolytope.Rd deleted file mode 100644 index e330af8b..00000000 --- a/man/Vpolytope.Rd +++ /dev/null @@ -1,17 +0,0 @@ -\name{Vpolytope} -\alias{Vpolytope} -\title{An \code{R} class to represent V-polytopes.} - -\description{ - A V-polytope is defined as the convex hull of \eqn{m} \eqn{d}-dimensional points which corresponds to its vertices. -} -\section{Fields}{ -\itemize{ - \item{\code{V} }{ A \eqn{m\times d} numerical matrix that contains the vertices row-wise} - - \item{\code{type} }{ An integer that declares the representation of the polytope. For V-representation the default value is 2.} - - \item{\code{dimension} }{ The dimension of the polytope.} - - \item{\code{volume} }{ The volume of the polytope, if it is known.} -}} diff --git a/man/VpolytopeIntersection-class.Rd b/man/VpolytopeIntersection-class.Rd new file mode 100644 index 00000000..ae6e9364 --- /dev/null +++ b/man/VpolytopeIntersection-class.Rd @@ -0,0 +1,27 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/VpolytopeIntersectionClass.R +\docType{class} +\name{VpolytopeIntersection-class} +\alias{VpolytopeIntersection-class} +\alias{VpolytopeIntersection} +\title{An R class to represent the intersection of two V-polytopes} +\description{ +An intersection of two V-polytopes is defined by the intersection of the two coresponding convex hulls. +} +\details{ +\describe{ + \item{V1}{An \eqn{m\times d} numerical matrix that contains the vertices of the first V-polytope (row-wise).} + + \item{V2}{An \eqn{q\times d} numerical matrix that contains the vertices of the second V-polytope (row-wise).} + + \item{volume}{The volume of the polytope if it is known, \eqn{NaN} otherwise by default.} + + \item{type}{A character with default value 'VpolytopeIntersection', to declare the representation of the polytope.} +} +} +\examples{ +P1 = gen_simplex(2,'V') +P2 = gen_cross(2,'V') +P = VpolytopeIntersection(V1 = P1@V, V2 = P2@V) + +} diff --git a/man/VpolytopeIntersection.Rd b/man/VpolytopeIntersection.Rd deleted file mode 100644 index f7a75c66..00000000 --- a/man/VpolytopeIntersection.Rd +++ /dev/null @@ -1,19 +0,0 @@ -\name{VpolytopeIntersection} -\alias{VpolytopeIntersection} -\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 corresponding convex hulls. -} -\section{Fields}{ -\itemize{ -\item{\code{V1} }{ The numerical matrix that contains the vertices of \eqn{P_1} row-wise.} - -\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 this polytope the default value is 4.} - -\item{\code{dimension} }{ The dimension of the polytope.} - -\item{\code{volume} }{ The volume of the polytope, if it is known.} -}} diff --git a/man/Zonotope-class.Rd b/man/Zonotope-class.Rd new file mode 100644 index 00000000..0e857d50 --- /dev/null +++ b/man/Zonotope-class.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ZonotopeClass.R +\docType{class} +\name{Zonotope-class} +\alias{Zonotope-class} +\alias{Zonotope} +\title{An R class to represent a Zonotope} +\description{ +A zonotope is a convex polytope defined by the Minkowski sum of \eqn{m} \eqn{d}-dimensional segments. +} +\details{ +\describe{ + \item{G}{An \eqn{m\times d} numerical matrix that contains the segments (or generators) row-wise} + + \item{volume}{The volume of the polytope if it is known, \eqn{NaN} otherwise by default.} + + \item{type}{A character with default value 'Zonotope', to declare the representation of the polytope.} +} +} +\examples{ +G = matrix(c(2,3,-1,7,0,0),ncol = 2, nrow = 3, byrow = TRUE) +P = Zonotope(G = G) + +} diff --git a/man/Zonotope.Rd b/man/Zonotope.Rd deleted file mode 100644 index 959fbed5..00000000 --- a/man/Zonotope.Rd +++ /dev/null @@ -1,17 +0,0 @@ -\name{Zonotope} -\alias{Zonotope} -\title{An \code{R} class to represent zonotopes.} - -\description{ -A zonotope is a convex polytope defined by the Minkowski sum of \eqn{m} \eqn{d}-dimensional segments. -} -\section{Fields}{ -\itemize{ -\item{\code{G} }{ A \eqn{m\times d} numerical matrix that contains the segments (or generators) row-wise} - -\item{\code{type} }{ An integer that declares the representation of the polytope. For zonotopes the default value is 3.} - -\item{\code{dimension} }{ The dimension of the polytope.} - -\item{\code{volume} }{ The volume of the polytope, if it is known.} -}} diff --git a/man/compute_indicators.Rd b/man/compute_indicators.Rd index c0eaadd9..d1d870a3 100644 --- a/man/compute_indicators.Rd +++ b/man/compute_indicators.Rd @@ -6,41 +6,35 @@ \usage{ compute_indicators( returns, - win_length = NULL, - m = NULL, - n = NULL, - nwarning = NULL, - ncrisis = NULL, - seed = NULL + parameters = list(win_length = 60, m = 100, n = 5e+05, nwarning = 60, ncrisis = 100) ) } \arguments{ \item{returns}{A \eqn{d}-dimensional vector that describes the direction of the first family of parallel hyperplanes.} -\item{win_length}{Optional. The length of the sliding window. The default value is 60.} - -\item{m}{Optional. The number of slices for the copula. The default value is 100.} - -\item{n}{Optional. The number of points to sample. The default value is \eqn{5\cdot 10^5}.} - -\item{nwarning}{Optional. The number of consecutive indicators larger than 1 required to declare a warning period. The default value is 60.} - -\item{ncrisis}{Optional. The number of consecutive indicators larger than 1 required to declare a crisis period. The default value is 100.} - -\item{seed}{Optional. A fixed seed for the number generator.} +\item{parameters}{A list to set a parameterization. +\itemize{ +\item{win_length }{ The length of the sliding window. The default value is 60.} +\item{m } { The number of slices for the copula. The default value is 100.} +\item{n }{ The number of points to sample. The default value is \eqn{5\cdot 10^5}.} +\item{nwarning }{ The number of consecutive indicators larger than 1 required to declare a warning period. The default value is 60.} +\item{ncrisis }{ The number of consecutive indicators larger than 1 required to declare a crisis period. The default value is 100.} +\item{seed }{ A fixed seed for the number generator.} +}} } \value{ A list that contains the indicators and the corresponding vector that label each time period with respect to the market state: a) normal, b) crisis, c) warning. } \description{ -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. +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) +market_states_and_indicators = compute_indicators(asset_returns, + parameters = list("win_length" = 10, "m" = 10, "n" = 10000, "nwarning" = 2, "ncrisis" = 3)) } \references{ diff --git a/man/exact_vol.Rd b/man/exact_vol.Rd index 59324740..33b58068 100644 --- a/man/exact_vol.Rd +++ b/man/exact_vol.Rd @@ -23,8 +23,8 @@ 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$new(V) +A = matrix(c(2,3,-1,7,0,0),ncol = 2, nrow = 3, byrow = TRUE) +P = Vpolytope(V = A) vol = exact_vol(P) } diff --git a/man/file_to_polytope.Rd b/man/file_to_polytope.Rd deleted file mode 100644 index 8ec843d5..00000000 --- a/man/file_to_polytope.Rd +++ /dev/null @@ -1,20 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/file_to_polytope.R -\name{file_to_polytope} -\alias{file_to_polytope} -\title{function to get an ine or an ext file and returns the corresponding polytope} -\usage{ -file_to_polytope(path, zonotope = FALSE) -} -\arguments{ -\item{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.} - -\item{zonotope}{A boolean parameter. It has to be TRUE when the path leads to an .ext file that describes a zonotope.} -} -\value{ -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 more details on those file formats see \url{https://github.com/GeomScale/volume_approximation/blob/develop/doc/cpp_interface.md#polytope-input}. -} diff --git a/man/gen_cross.Rd b/man/gen_cross.Rd index 53e84c78..ea041442 100644 --- a/man/gen_cross.Rd +++ b/man/gen_cross.Rd @@ -4,12 +4,12 @@ \alias{gen_cross} \title{Generator function for cross polytopes} \usage{ -gen_cross(dimension, representation) +gen_cross(dimension, representation = "H") } \arguments{ \item{dimension}{The dimension of the cross polytope.} -\item{representation}{A string to declare the representation. It has to be \code{'H'} for H-representation or \code{'V'} for V-representation.} +\item{representation}{A string to declare the representation. It has to be \code{'H'} for H-representation or \code{'V'} for V-representation. Default valus is 'H'.} } \value{ A polytope class representing a cross polytope in H- or V-representation. diff --git a/man/gen_cube.Rd b/man/gen_cube.Rd index efb486d6..dc97ae5a 100644 --- a/man/gen_cube.Rd +++ b/man/gen_cube.Rd @@ -4,12 +4,12 @@ \alias{gen_cube} \title{Generator function for hypercubes} \usage{ -gen_cube(dimension, representation) +gen_cube(dimension, representation = "H") } \arguments{ \item{dimension}{The dimension of the hypercube} -\item{representation}{A string to declare the representation. It has to be \code{'H'} for H-representation or \code{'V'} for V-representation.} +\item{representation}{A string to declare the representation. It has to be \code{'H'} for H-representation or \code{'V'} for V-representation. Default valus is 'H'.} } \value{ A polytope class representing the unit \eqn{d}-dimensional hypercube in H- or V-representation. diff --git a/man/gen_rand_hpoly.Rd b/man/gen_rand_hpoly.Rd index e4db1c9e..85dd1757 100644 --- a/man/gen_rand_hpoly.Rd +++ b/man/gen_rand_hpoly.Rd @@ -4,14 +4,18 @@ \alias{gen_rand_hpoly} \title{Generator function for random H-polytopes} \usage{ -gen_rand_hpoly(dimension, nfacets, seed = NULL) +gen_rand_hpoly(dimension, nfacets, generator = list(constants = "sphere")) } \arguments{ \item{dimension}{The dimension of the convex polytope.} \item{nfacets}{The number of the facets.} -\item{seed}{Optional. A fixed seed for the generator.} +\item{generator}{A list that could contain two elements. +\itemize{ +\item{constants }{ To declare how to set the constants \eqn{b_i} for each facets: (i) 'sphere', each hyperplane is tangent to the hypersphere of radius 10, (ii) 'uniform' for each \eqn{b_i} the generator picks a uniform number from \eqn{(0,1)}. The defalut value is 'sphere'.} +\item{seed }{ Optional. A fixed seed for the number generator.} +}} } \value{ A polytope class representing a H-polytope. diff --git a/man/gen_rand_vpoly.Rd b/man/gen_rand_vpoly.Rd index e2c3c097..e34125a5 100644 --- a/man/gen_rand_vpoly.Rd +++ b/man/gen_rand_vpoly.Rd @@ -4,16 +4,18 @@ \alias{gen_rand_vpoly} \title{Generator function for random V-polytopes} \usage{ -gen_rand_vpoly(dimension, nvertices, generator = NULL, seed = NULL) +gen_rand_vpoly(dimension, nvertices, generator = list(body = "sphere")) } \arguments{ \item{dimension}{The dimension of the convex polytope.} \item{nvertices}{The number of the vertices.} -\item{generator}{The body that the generator samples uniformly the vertices from: (a) 'cube' or (b) 'sphere'.} - -\item{seed}{Optional. A fixed seed for the generator.} +\item{generator}{A list that could contain two elements. +\itemize{ +\item{body }{ the body that the generator samples uniformly the vertices from: (i) 'cube' or (ii) 'sphere', the default value is 'sphere'.} +\item{seed }{ Optional. A fixed seed for the number generator.} +}} } \value{ A polytope class representing a V-polytope. diff --git a/man/gen_rand_zonotope.Rd b/man/gen_rand_zonotope.Rd index ab971d35..41497ef2 100644 --- a/man/gen_rand_zonotope.Rd +++ b/man/gen_rand_zonotope.Rd @@ -4,16 +4,22 @@ \alias{gen_rand_zonotope} \title{Generator function for zonotopes} \usage{ -gen_rand_zonotope(dimension, nsegments, generator = NULL, seed = NULL) +gen_rand_zonotope( + dimension, + nsegments, + generator = list(distribution = "uniform") +) } \arguments{ \item{dimension}{The dimension of the zonotope.} \item{nsegments}{The number of segments that generate the zonotope.} -\item{generator}{The distribution to pick the length of each segment from \eqn{[0,100]}: (a) 'uniform', (b) 'gaussian' or (c) 'exponential'.} - -\item{seed}{Optional. A fixed seed for the generator.} +\item{generator}{A list that could contain two elements. +\itemize{ +\item{distribution }{ the distribution to pick the length of each segment from \eqn{[0,100]}: (i) 'uniform', (ii) 'gaussian' or (iii) 'exponential', the default value is 'uniform.} +\item {seed }{ Optional. A fixed seed for the number generator.} +}} } \value{ A polytope class representing a zonotope. diff --git a/man/gen_simplex.Rd b/man/gen_simplex.Rd index bed04a3a..0bbcb155 100644 --- a/man/gen_simplex.Rd +++ b/man/gen_simplex.Rd @@ -4,12 +4,12 @@ \alias{gen_simplex} \title{Generator function for simplices} \usage{ -gen_simplex(dimension, representation) +gen_simplex(dimension, representation = "H") } \arguments{ \item{dimension}{The dimension of the unit simplex.} -\item{representation}{A string to declare the representation. It has to be \code{'H'} for H-representation or \code{'V'} for V-representation.} +\item{representation}{A string to declare the representation. It has to be \code{'H'} for H-representation or \code{'V'} for V-representation. Default valus is 'H'.} } \value{ A polytope class representing the \eqn{d}-dimensional unit simplex in H- or V-representation. diff --git a/man/load_sdpa_format_file.Rd b/man/load_sdpa_format_file.Rd new file mode 100644 index 00000000..349c0dfc --- /dev/null +++ b/man/load_sdpa_format_file.Rd @@ -0,0 +1,18 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{load_sdpa_format_file} +\alias{load_sdpa_format_file} +\title{An internal Rccp function to read a SDPA format file} +\usage{ +load_sdpa_format_file(input_file = NULL) +} +\arguments{ +\item{input_file}{Name of the input file} +} +\value{ +A list with two named items: an item "matrices" which is a list of the matrices and an vector "objFunction" +} +\description{ +An internal Rccp function to read a SDPA format file +} +\keyword{internal} diff --git a/man/readSdpaFormatFile.Rd b/man/read_sdpa_format_file.Rd similarity index 81% rename from man/readSdpaFormatFile.Rd rename to man/read_sdpa_format_file.Rd index 878a1196..50ca8ae2 100644 --- a/man/readSdpaFormatFile.Rd +++ b/man/read_sdpa_format_file.Rd @@ -1,10 +1,10 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/read_sdpa_file.R -\name{readSdpaFormatFile} -\alias{readSdpaFormatFile} +\name{read_sdpa_format_file} +\alias{read_sdpa_format_file} \title{Read a SDPA format file} \usage{ -readSdpaFormatFile(path) +read_sdpa_format_file(path) } \arguments{ \item{path}{Name of the input file} @@ -18,7 +18,7 @@ the linear matrix inequality in the input file, and the objective function. } \examples{ path = system.file('extdata', package = 'volesti') -l = readSdpaFormatFile(paste0(path,'/sdpa_n2m3.txt')) +l = read_sdpa_format_file(paste0(path,'/sdpa_n2m3.txt')) Spectrahedron = l$spectrahedron objFunction = l$objFunction } diff --git a/man/rotate_polytope.Rd b/man/rotate_polytope.Rd index 70c9c3ad..1368e515 100644 --- a/man/rotate_polytope.Rd +++ b/man/rotate_polytope.Rd @@ -4,14 +4,12 @@ \alias{rotate_polytope} \title{Apply a random rotation to a convex polytope (H-polytope, V-polytope, zonotope or intersection of two V-polytopes)} \usage{ -rotate_polytope(P, T = NULL, seed = NULL) +rotate_polytope(P, rotation = list()) } \arguments{ \item{P}{A convex polytope. It is an object from class (a) Hpolytope, (b) Vpolytope, (c) Zonotope, (d) intersection of two V-polytopes.} -\item{T}{Optional. A \eqn{d\times d} rotation matrix.} - -\item{seed}{Optional. A fixed seed for the random linear map generator.} +\item{rotation}{A list that contains (a) the rotation matrix T and (b) the 'seed' to set a spesific seed for the number generator.} } \value{ A list that contains the rotated polytope and the matrix \eqn{T} of the linear transformation. @@ -20,7 +18,7 @@ A list that contains the rotated polytope and the matrix \eqn{T} of the linear t Given a convex H- or V- polytope or a zonotope or an intersection of two V-polytopes as input, this function applies (a) a random rotation or (b) a given rotation by an input matrix \eqn{T}. } \details{ -Let \eqn{P} be the given polytope and \eqn{Q} the rotated one and \eqn{T} be the matrix of the linear transformation. +Let \eqn{P} be the given polytope and \eqn{Q} the rotated one and \eqn{T} be the matrix of the linear transformation. \itemize{ \item{If \eqn{P} is in H-representation and \eqn{A} is the matrix that contains the normal vectors of the facets of \eqn{Q} then \eqn{AT} contains the normal vactors of the facets of \eqn{P}.} \item{If \eqn{P} is in V-representation and \eqn{V} is the matrix that contains column-wise the vertices of \eqn{Q} then \eqn{T^TV} contains the vertices of \eqn{P}.} @@ -30,7 +28,7 @@ Let \eqn{P} be the given polytope and \eqn{Q} the rotated one and \eqn{T} be the } \examples{ # rotate a H-polytope (2d unit simplex) -P = gen_simplex(2,'H') +P = gen_simplex(2, 'H') poly_matrix_list = rotate_polytope(P) # rotate a V-polytope (3d cube) @@ -38,6 +36,6 @@ P = gen_cube(3, 'V') poly_matrix_list = rotate_polytope(P) # rotate a 5-dimensional zonotope defined by the Minkowski sum of 15 segments -Z = gen_rand_zonotope(3,6) +Z = gen_rand_zonotope(3, 6) poly_matrix_list = rotate_polytope(Z) } diff --git a/man/round_polytope.Rd b/man/round_polytope.Rd index 787e40fb..04b91bf3 100644 --- a/man/round_polytope.Rd +++ b/man/round_polytope.Rd @@ -4,14 +4,16 @@ \alias{round_polytope} \title{Apply rounding to a convex polytope (H-polytope, V-polytope or a zonotope)} \usage{ -round_polytope(P, method = NULL, seed = NULL) +round_polytope(P, settings = list()) } \arguments{ \item{P}{A convex polytope. It is an object from class (a) Hpolytope or (b) Vpolytope or (c) Zonotope.} -\item{method}{Optional. The method to use for rounding, a) \code{'min_ellipsoid'} for the method based on mimimmum volume enclosing ellipsoid of a dataset, b) \code{'max_ellipsoid'} for the method based on maximum volume enclosed ellipsoid, (c) \code{'isotropy'} for the method based on svd decomposition. The default method is \code{'mee'} for all the representations.} - -\item{seed}{Optional. A fixed seed for the number generator.} +\item{settings}{Optional. A list of settings. +\itemize{ +\item{\code{method} }{ The method to use for rounding, a) \code{'min_ellipsoid'} for the method based on mimimmum volume enclosing ellipsoid of a dataset, b) \code{'max_ellipsoid'} for the method based on maximum volume enclosed ellipsoid, (c) \code{'isotropy'} for the method based on svd decomposition. The default method is \code{'mee'} for all the representations.} +\item{\code{seed} }{ Optional. A fixed seed for the number generator.} +}} } \value{ 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. diff --git a/man/sample_points.Rd b/man/sample_points.Rd index f5ea6f2b..ed0cc961 100644 --- a/man/sample_points.Rd +++ b/man/sample_points.Rd @@ -13,25 +13,26 @@ sample_points(P, n, random_walk = NULL, distribution = NULL, seed = NULL) \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{'dikin'} for dikin walk, vi) \code{'vaidya'} for vaidya walk, vii) \code{'john'} for john walk, viii) \code{'BCDHR'} boundary sampling by keeping the extreme points of CDHR or ix) \code{'BRDHR'} boundary sampling by keeping the extreme points of RDHR x) \code{'HMC'} for Hamiltonian Monte Carlo (logconcave) xi) \code{'ULD'} for Underdamped Langevin Dynamics using the Randomized Midpoint Method. The default walk is \code{'aBiW'} for the uniform distribution or \code{'CDHR'} for the Gaussian distribution and H-polytopes and \code{'BiW'} or \code{'RDHR'} for the same distributions and V-polytopes and zonotopes.} +\item{\code{walk} }{ A string to declare the random walk: i) \code{'CDHR'} for Coordinate Directions Hit-and-Run, ii) \code{'RDHR'} for Random Directions Hit-and-Run, iii) \code{'BaW'} for Ball Walk, iv) \code{'BiW'} for Billiard walk, v) \code{'dikin'} for dikin walk, vi) \code{'vaidya'} for vaidya walk, vii) \code{'john'} for john walk, viii) \code{'BCDHR'} boundary sampling by keeping the extreme points of CDHR or ix) \code{'BRDHR'} boundary sampling by keeping the extreme points of RDHR x) \code{'HMC'} for Hamiltonian Monte Carlo (logconcave densities) xi) \code{'ULD'} for Underdamped Langevin Dynamics using the Randomized Midpoint Method xii) \code{'ExactHMC'} for exact Hamiltonian Monte Carlo with reflections (spherical Gaussian or exponential distribution). The default walk is \code{'aBiW'} for the uniform distribution or \code{'CDHR'} for the Gaussian distribution and H-polytopes and \code{'BiW'} or \code{'RDHR'} for the same distributions and V-polytopes and zonotopes.} \item{\code{walk_length} }{ The number of the steps per generated point for the random walk. The default value is \eqn{1}.} \item{\code{nburns} }{ The number of points to burn before start sampling. The default value is \eqn{1}.} \item{\code{starting_point} }{ A \eqn{d}-dimensional numerical vector that declares a starting point in the interior of the polytope for the random walk. The default choice is the center of the ball as that one computed by the function \code{inner_ball()}.} \item{\code{BaW_rad} }{ The radius for the ball walk.} \item{\code{L} }{ The maximum length of the billiard trajectory or the radius for the step of dikin, vaidya or john walk.} -\item{\code{solver}} {Specify ODE solver for logconcave sampling. Options are i) leapfrog, ii) euler iii) runge-kutta iv) richardson} -\item{\code{step_size} {Optionally chosen step size for logconcave sampling. Defaults to a theoretical value if not provided.}} +\item{\code{solver} }{ Specify ODE solver for logconcave sampling. Options are i) leapfrog, ii) euler iii) runge-kutta iv) richardson} +\item{\code{step_size} }{ Optionally chosen step size for logconcave sampling. Defaults to a theoretical value if not provided.} }} \item{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. c) Logconcave with form proportional to exp(-f(x)) where f(x) is L-smooth and m-strongly-convex. } -\item{\code{variance} }{ The variance of the multidimensional spherical gaussian. The default value is 1.} +\item{\code{density} }{ A string: (a) \code{'uniform'} for the uniform distribution or b) \code{'gaussian'} for the multidimensional spherical distribution c) \code{logconcave} with form proportional to exp(-f(x)) where f(x) is L-smooth and m-strongly-convex d) \code{'exponential'} for the exponential distribution. The default target distribution is the uniform distribution.} +\item{\code{variance} }{ The variance of the multidimensional spherical gaussian or the exponential distribution. 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()}.} -\item{\code{L_}} { Smoothness constant (for logconcave). } -\item{\code{m}} { Strong-convexity constant (for logconcave). } -\item{\code{negative_logprob}} { Negative log-probability (for logconcave). } -\item{\code{negative_logprob_gradient}} { Negative log-probability gradient (for logconcave). } +\item{\code{bias} }{ The bias vector for the exponential distribution. The default vector is \eqn{c_1 = 1} and \eqn{c_i = 0} for \eqn{i \neq 1}.} +\item{\code{L_} }{ Smoothness constant (for logconcave). } +\item{\code{m} }{ Strong-convexity constant (for logconcave). } +\item{\code{negative_logprob} }{ Negative log-probability (for logconcave). } +\item{\code{negative_logprob_gradient} }{ Negative log-probability gradient (for logconcave). } }} \item{seed}{Optional. A fixed seed for the number generator.} @@ -50,7 +51,7 @@ points = sample_points(P, n = 100, random_walk = list("walk" = "BaW", "walk_leng # gaussian distribution from the 2d unit simplex in H-representation with variance = 2 A = matrix(c(-1,0,0,-1,1,1), ncol=2, nrow=3, byrow=TRUE) b = c(0,0,1) -P = Hpolytope$new(A,b) +P = Hpolytope(A = A, b = b) points = sample_points(P, n = 100, distribution = list("density" = "gaussian", "variance" = 2)) # uniform points from the boundary of a 2-dimensional random H-polytope @@ -75,4 +76,7 @@ points = sample_points(P, n = 100, random_walk = list("walk" = "BRDHR")) \cite{Shen, Ruoqi, and Yin Tat Lee, \dQuote{"The randomized midpoint method for log-concave sampling.",} \emph{Advances in Neural Information Processing Systems}, 2019.} + +\cite{Augustin Chevallier, Sylvain Pion, Frederic Cazals, +\dQuote{"Hamiltonian Monte Carlo with boundary reflections, and application to polytope volume calculations,"} \emph{Research Report preprint hal-01919855}, 2018.} } diff --git a/man/volume.Rd b/man/volume.Rd index cd325c87..d00d797d 100644 --- a/man/volume.Rd +++ b/man/volume.Rd @@ -4,7 +4,7 @@ \alias{volume} \title{The main function for volume approximation of a convex Polytope (H-polytope, V-polytope, zonotope or intersection of two V-polytopes). It returns a list with two elements: (a) the logarithm of the estimated volume and (b) the estimated volume} \usage{ -volume(P, settings = NULL, rounding = NULL, seed = NULL) +volume(P, settings = NULL, rounding = NULL) } \arguments{ \item{P}{A convex polytope. It is an object from class a) Hpolytope or b) Vpolytope or c) Zonotope or d) VpolytopeIntersection.} @@ -17,6 +17,7 @@ volume(P, settings = NULL, rounding = NULL, seed = NULL) \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{250} for CB with BiW and \eqn{400+3d^2} for CB and any other random walk and \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.} }} \item{rounding}{Optional. A string parameter to request a rounding method to be applied in the input polytope before volume computation: a) \code{'min_ellipsoid'}, b) \code{'svd'}, c) \code{'max_ellipsoid'} and d) \code{'none'} for no rounding.} diff --git a/man/writeSdpaFormatFile.Rd b/man/writeSdpaFormatFile.Rd index 7c791dce..f70cac91 100644 --- a/man/writeSdpaFormatFile.Rd +++ b/man/writeSdpaFormatFile.Rd @@ -27,7 +27,7 @@ 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$new(lmi); +S = Spectrahedron(matrices = lmi) objFunction = c(1,1) writeSdpaFormatFile(S, objFunction, "output.txt") } diff --git a/man/write_sdpa_format_file.Rd b/man/write_sdpa_format_file.Rd new file mode 100644 index 00000000..23ec6714 --- /dev/null +++ b/man/write_sdpa_format_file.Rd @@ -0,0 +1,30 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{write_sdpa_format_file} +\alias{write_sdpa_format_file} +\title{Write a SDPA format file} +\usage{ +write_sdpa_format_file(spectrahedron, objective_function, output_file) +} +\arguments{ +\item{spectrahedron}{A spectrahedron in n dimensions; must be an object of class Spectrahedron} + +\item{objective_function}{A numerical vector of length n} + +\item{output_file}{Name of the output file} +} +\description{ +Outputs a spectrahedron (the matrices defining a linear matrix inequality) and a vector (the objective function) +to a SDPA format 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") +} +} diff --git a/man/zonotope_approximation.Rd b/man/zonotope_approximation.Rd index 152a0b0a..b06743a6 100644 --- a/man/zonotope_approximation.Rd +++ b/man/zonotope_approximation.Rd @@ -4,7 +4,11 @@ \alias{zonotope_approximation} \title{A function to over-approximate a zonotope with PCA method and to evaluate the approximation by computing a ratio of fitness.} \usage{ -zonotope_approximation(Z, fit_ratio = NULL, settings = NULL, seed = NULL) +zonotope_approximation( + Z, + fit_ratio = FALSE, + settings = list(error = 0.1, walk_length = 1, win_len = 250, hpoly = FALSE) +) } \arguments{ \item{Z}{A zonotope.} @@ -15,11 +19,10 @@ zonotope_approximation(Z, fit_ratio = NULL, settings = NULL, seed = NULL) \itemize{ \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{win_len} }{ The length of the sliding window for CB algorithm. The default value is \eqn{250}.} \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{\code{seed} }{ Optional. A fixed seed for the number generator.} }} - -\item{seed}{Optional. A fixed seed for the number generator.} } \value{ A list that contains the approximation body in H-representation and the ratio of fitness @@ -29,7 +32,7 @@ For the evaluation of the PCA method the exact volume of the approximation body } \examples{ # over-approximate a 2-dimensional zonotope with 10 generators and compute the ratio of fitness -Z = gen_rand_zonotope(2,12) +Z = gen_rand_zonotope(2, 10) retList = zonotope_approximation(Z = Z) } diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index bc377c76..93e873c4 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -52,12 +52,12 @@ BEGIN_RCPP END_RCPP } // exact_vol -double exact_vol(Rcpp::Nullable P); +double exact_vol(Rcpp::Reference P); RcppExport SEXP _volesti_exact_vol(SEXP PSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< Rcpp::Nullable >::type P(PSEXP); + Rcpp::traits::input_parameter< Rcpp::Reference >::type P(PSEXP); rcpp_result_gen = Rcpp::wrap(exact_vol(P)); return rcpp_result_gen; END_RCPP @@ -99,6 +99,17 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// load_sdpa_format_file +Rcpp::List load_sdpa_format_file(Rcpp::Nullable input_file); +RcppExport SEXP _volesti_load_sdpa_format_file(SEXP input_fileSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< Rcpp::Nullable >::type input_file(input_fileSEXP); + rcpp_result_gen = Rcpp::wrap(load_sdpa_format_file(input_file)); + return rcpp_result_gen; +END_RCPP +} // ode_solve Rcpp::List ode_solve(Rcpp::Nullable n, Rcpp::Nullable step_size, Rcpp::Nullable order, Rcpp::Nullable dimension, Rcpp::Nullable initial_time, Rcpp::Function F, Rcpp::Nullable method, Rcpp::Nullable domains, Rcpp::Nullable initial_conditions); RcppExport SEXP _volesti_ode_solve(SEXP nSEXP, SEXP step_sizeSEXP, SEXP orderSEXP, SEXP dimensionSEXP, SEXP initial_timeSEXP, SEXP FSEXP, SEXP methodSEXP, SEXP domainsSEXP, SEXP initial_conditionsSEXP) { @@ -198,12 +209,12 @@ BEGIN_RCPP END_RCPP } // sample_points -Rcpp::NumericMatrix sample_points(Rcpp::Nullable P, Rcpp::Nullable n, Rcpp::Nullable random_walk, Rcpp::Nullable distribution, Rcpp::Nullable seed); +Rcpp::NumericMatrix sample_points(Rcpp::Reference P, Rcpp::Nullable n, Rcpp::Nullable random_walk, Rcpp::Nullable distribution, Rcpp::Nullable seed); RcppExport SEXP _volesti_sample_points(SEXP PSEXP, SEXP nSEXP, SEXP random_walkSEXP, SEXP distributionSEXP, SEXP seedSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< Rcpp::Nullable >::type P(PSEXP); + Rcpp::traits::input_parameter< Rcpp::Reference >::type P(PSEXP); Rcpp::traits::input_parameter< Rcpp::Nullable >::type n(nSEXP); Rcpp::traits::input_parameter< Rcpp::Nullable >::type random_walk(random_walkSEXP); Rcpp::traits::input_parameter< Rcpp::Nullable >::type distribution(distributionSEXP); @@ -213,11 +224,11 @@ BEGIN_RCPP END_RCPP } // writeSdpaFormatFile -void writeSdpaFormatFile(Rcpp::Nullable spectrahedron, Rcpp::Nullable objectiveFunction, Rcpp::Nullable outputFile); +void writeSdpaFormatFile(Rcpp::Reference spectrahedron, Rcpp::Nullable objectiveFunction, Rcpp::Nullable outputFile); RcppExport SEXP _volesti_writeSdpaFormatFile(SEXP spectrahedronSEXP, SEXP objectiveFunctionSEXP, SEXP outputFileSEXP) { BEGIN_RCPP Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< Rcpp::Nullable >::type spectrahedron(spectrahedronSEXP); + Rcpp::traits::input_parameter< Rcpp::Reference >::type spectrahedron(spectrahedronSEXP); Rcpp::traits::input_parameter< Rcpp::Nullable >::type objectiveFunction(objectiveFunctionSEXP); Rcpp::traits::input_parameter< Rcpp::Nullable >::type outputFile(outputFileSEXP); writeSdpaFormatFile(spectrahedron, objectiveFunction, outputFile); @@ -236,19 +247,30 @@ BEGIN_RCPP END_RCPP } // volume -Rcpp::List volume(Rcpp::Reference P, Rcpp::Nullable settings, Rcpp::Nullable rounding, Rcpp::Nullable seed); -RcppExport SEXP _volesti_volume(SEXP PSEXP, SEXP settingsSEXP, SEXP roundingSEXP, SEXP seedSEXP) { +Rcpp::List volume(Rcpp::Reference P, Rcpp::Nullable settings, Rcpp::Nullable rounding); +RcppExport SEXP _volesti_volume(SEXP PSEXP, SEXP settingsSEXP, SEXP roundingSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; Rcpp::traits::input_parameter< Rcpp::Reference >::type P(PSEXP); Rcpp::traits::input_parameter< Rcpp::Nullable >::type settings(settingsSEXP); Rcpp::traits::input_parameter< Rcpp::Nullable >::type rounding(roundingSEXP); - Rcpp::traits::input_parameter< Rcpp::Nullable >::type seed(seedSEXP); - rcpp_result_gen = Rcpp::wrap(volume(P, settings, rounding, seed)); + rcpp_result_gen = Rcpp::wrap(volume(P, settings, rounding)); return rcpp_result_gen; END_RCPP } +// write_sdpa_format_file +void write_sdpa_format_file(Rcpp::Reference spectrahedron, Rcpp::NumericVector objective_function, std::string output_file); +RcppExport SEXP _volesti_write_sdpa_format_file(SEXP spectrahedronSEXP, SEXP objective_functionSEXP, SEXP output_fileSEXP) { +BEGIN_RCPP + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< Rcpp::Reference >::type spectrahedron(spectrahedronSEXP); + Rcpp::traits::input_parameter< Rcpp::NumericVector >::type objective_function(objective_functionSEXP); + Rcpp::traits::input_parameter< std::string >::type output_file(output_fileSEXP); + write_sdpa_format_file(spectrahedron, objective_function, output_file); + return R_NilValue; +END_RCPP +} // zono_approx Rcpp::List zono_approx(Rcpp::Reference Z, Rcpp::Nullable fit_ratio, Rcpp::Nullable settings, Rcpp::Nullable seed); RcppExport SEXP _volesti_zono_approx(SEXP ZSEXP, SEXP fit_ratioSEXP, SEXP settingsSEXP, SEXP seedSEXP) { @@ -272,6 +294,7 @@ static const R_CallMethodDef CallEntries[] = { {"_volesti_frustum_of_simplex", (DL_FUNC) &_volesti_frustum_of_simplex, 2}, {"_volesti_geweke", (DL_FUNC) &_volesti_geweke, 3}, {"_volesti_inner_ball", (DL_FUNC) &_volesti_inner_ball, 2}, + {"_volesti_load_sdpa_format_file", (DL_FUNC) &_volesti_load_sdpa_format_file, 1}, {"_volesti_ode_solve", (DL_FUNC) &_volesti_ode_solve, 9}, {"_volesti_poly_gen", (DL_FUNC) &_volesti_poly_gen, 6}, {"_volesti_psrf_multivariate", (DL_FUNC) &_volesti_psrf_multivariate, 1}, @@ -282,7 +305,8 @@ static const R_CallMethodDef CallEntries[] = { {"_volesti_sample_points", (DL_FUNC) &_volesti_sample_points, 5}, {"_volesti_writeSdpaFormatFile", (DL_FUNC) &_volesti_writeSdpaFormatFile, 3}, {"_volesti_loadSdpaFormatFile", (DL_FUNC) &_volesti_loadSdpaFormatFile, 1}, - {"_volesti_volume", (DL_FUNC) &_volesti_volume, 4}, + {"_volesti_volume", (DL_FUNC) &_volesti_volume, 3}, + {"_volesti_write_sdpa_format_file", (DL_FUNC) &_volesti_write_sdpa_format_file, 3}, {"_volesti_zono_approx", (DL_FUNC) &_volesti_zono_approx, 4}, {NULL, NULL, 0} }; diff --git a/src/copula.cpp b/src/copula.cpp index 0e8069ea..8b9fd586 100644 --- a/src/copula.cpp +++ b/src/copula.cpp @@ -7,7 +7,6 @@ //Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 program. - #include #include #include diff --git a/src/exact_vol.cpp b/src/exact_vol.cpp index d60d8fa4..89bc3a36 100644 --- a/src/exact_vol.cpp +++ b/src/exact_vol.cpp @@ -43,8 +43,8 @@ FT factorial(FT n) //' 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$new(V) +//' A = matrix(c(2,3,-1,7,0,0),ncol = 2, nrow = 3, byrow = TRUE) +//' P = Vpolytope(V = A) //' vol = exact_vol(P) //' } //' @@ -53,7 +53,7 @@ FT factorial(FT n) //' vol = exact_vol(P) //' @export // [[Rcpp::export]] -double exact_vol(Rcpp::Nullable P) { +double exact_vol(Rcpp::Reference P) { typedef double NT; typedef Cartesian Kernel; @@ -61,21 +61,38 @@ double exact_vol(Rcpp::Nullable P) { typedef Eigen::Matrix VT; typedef Eigen::Matrix MT; - if (NT(Rcpp::as(P).field("volume")) > 0.0) { - return NT(Rcpp::as(P).field("volume")); + if (NT(P.slot("volume")) > 0.0) { + return NT(P.slot("volume")); + } + + int type; + int dim; + std::string type_str = Rcpp::as(P.slot("type")); + + if (type_str.compare(std::string("Hpolytope")) == 0) { + dim = Rcpp::as(P.slot("A")).cols(); + type = 1; + } else if (type_str.compare(std::string("Vpolytope")) == 0) { + dim = Rcpp::as(P.slot("V")).cols(); + type = 2; + } else if (type_str.compare(std::string("Zonotope")) == 0) { + dim = Rcpp::as(P.slot("G")).cols(); + type = 3; + } else if (type_str.compare(std::string("VpolytopeIntersection")) == 0) { + dim = Rcpp::as(P.slot("V1")).cols(); + type = 4; + } else { + throw Rcpp::exception("Unknown polytope representation!"); } - int type = Rcpp::as(P).field("type"), dim; NT vol; if (type == 2) { - dim = Rcpp::as(P).field("dimension"); - - if (Rcpp::as(Rcpp::as(P).field("V")).rows() == - Rcpp::as(Rcpp::as(P).field("V")).cols() + 1) { + if (Rcpp::as(P.slot("V")).rows() == + Rcpp::as(P.slot("V")).cols() + 1) { - MT V = Rcpp::as(Rcpp::as(P).field("V")).transpose(), V2(dim,dim); + MT V = Rcpp::as(P.slot("V")).transpose(), V2(dim,dim); VT v0 = V.col(dim); for (int i = 0; i < dim; ++i) { @@ -90,10 +107,8 @@ double exact_vol(Rcpp::Nullable P) { } else if (type == 3) { typedef Zonotope zonotope; - dim = Rcpp::as(P).field("dimension"); - zonotope ZP(dim, Rcpp::as(Rcpp::as(P).field("G")), - VT::Ones(Rcpp::as(Rcpp::as(P).field("G")).rows())); + zonotope ZP(dim, Rcpp::as(P.slot("G")), VT::Ones(Rcpp::as(P.slot("G")).rows())); vol = exact_zonotope_vol(ZP); } else { diff --git a/src/inner_ball.cpp b/src/inner_ball.cpp index 8472ad3b..c69884e6 100644 --- a/src/inner_ball.cpp +++ b/src/inner_ball.cpp @@ -32,7 +32,7 @@ //' ball_vec = inner_ball(P, lpsolve = TRUE) //' @export // [[Rcpp::export]] -Rcpp::NumericVector inner_ball(Rcpp::Reference P, +Rcpp::NumericVector inner_ball(Rcpp::Reference P, Rcpp::Nullable lpsolve = R_NilValue) { typedef double NT; @@ -45,15 +45,34 @@ Rcpp::NumericVector inner_ball(Rcpp::Reference P, typedef IntersectionOfVpoly InterVP; typedef Eigen::Matrix VT; typedef Eigen::Matrix MT; - unsigned int n = P.field("dimension"), type = P.field("type"); std::pair InnerBall; bool lp_solve = (lpsolve.isNotNull()) ? Rcpp::as(lpsolve) : false; + unsigned int n; + unsigned int type; + std::string type_str = Rcpp::as(P.slot("type")); + + if (type_str.compare(std::string("Hpolytope")) == 0) { + n = Rcpp::as(P.slot("A")).cols(); + type = 1; + } else if (type_str.compare(std::string("Vpolytope")) == 0) { + n = Rcpp::as(P.slot("V")).cols(); + type = 2; + } else if (type_str.compare(std::string("Zonotope")) == 0) { + n = Rcpp::as(P.slot("G")).cols(); + type = 3; + } else if (type_str.compare(std::string("VpolytopeIntersection")) == 0) { + n = Rcpp::as(P.slot("V1")).cols(); + type = 4; + } else { + throw Rcpp::exception("Unknown polytope representation!"); + } + switch (type) { case 1: { // Hpolytope - Hpolytope HP(n, Rcpp::as(P.field("A")), Rcpp::as(P.field("b"))); + Hpolytope HP(n, Rcpp::as(P.slot("A")), Rcpp::as(P.slot("b"))); if (lp_solve) { InnerBall = ComputeChebychevBall(HP.get_mat(), HP.get_vec()); } else { @@ -64,22 +83,22 @@ Rcpp::NumericVector inner_ball(Rcpp::Reference P, } case 2: { // Vpolytope - Vpolytope VP(n, Rcpp::as(P.field("V")), VT::Ones(Rcpp::as(P.field("V")).rows())); + Vpolytope VP(n, Rcpp::as(P.slot("V")), VT::Ones(Rcpp::as(P.slot("V")).rows())); InnerBall = VP.ComputeInnerBall(); if (InnerBall.second < 0.0) throw Rcpp::exception("Unable to compute a feasible point."); break; } case 3: { // Zonotope - zonotope ZP(n, Rcpp::as(P.field("G")), VT::Ones(Rcpp::as(P.field("G")).rows())); + zonotope ZP(n, Rcpp::as(P.slot("G")), VT::Ones(Rcpp::as(P.slot("G")).rows())); InnerBall = ZP.ComputeInnerBall(); if (InnerBall.second < 0.0) throw Rcpp::exception("Unable to compute a feasible point."); break; } case 4: { // Intersection of two V-polytopes - Vpolytope VP1(n, Rcpp::as(P.field("V1")), VT::Ones(Rcpp::as(P.field("V1")).rows())); - Vpolytope VP2(n, Rcpp::as(P.field("V2")), VT::Ones(Rcpp::as(P.field("V2")).rows())); + Vpolytope VP1(n, Rcpp::as(P.slot("V1")), VT::Ones(Rcpp::as(P.slot("V1")).rows())); + Vpolytope VP2(n, Rcpp::as(P.slot("V2")), VT::Ones(Rcpp::as(P.slot("V2")).rows())); InterVP VPcVP(VP1, VP2); if (!VPcVP.is_feasible()) throw Rcpp::exception("Empty set!"); InnerBall = VPcVP.ComputeInnerBall(); diff --git a/src/load_sdpa_format_file.cpp b/src/load_sdpa_format_file.cpp new file mode 100644 index 00000000..bb1e00cd --- /dev/null +++ b/src/load_sdpa_format_file.cpp @@ -0,0 +1,62 @@ +// VolEsti (volume computation and sampling library) + +// Copyright (c) 2012-2020 Vissarion Fisikopoulos +// Copyright (c) 2020 Apostolos Chalkis + +//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. + +// Licensed under GNU LGPL.3, see LICENCE file + +#include +#include +#include +#include "cartesian_geom/cartesian_kernel.h" +#include +#include +#include +#include +#include "convex_bodies/spectrahedra/LMI.h" +#include "convex_bodies/spectrahedra/spectrahedron.h" +#include "SDPAFormatManager.h" + +//' 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" +// [[Rcpp::export]] +Rcpp::List load_sdpa_format_file(Rcpp::Nullable input_file = R_NilValue) { + + typedef double NT; + typedef Eigen::Matrix VT; + typedef Eigen::Matrix MT; + typedef Cartesian Kernel; + typedef typename Kernel::Point Point; + typedef boost::mt19937 RNGType; + typedef LMI LMI; + typedef Spectrahedron Spectrahedron; + + Spectrahedron _spectrahedron; + Point c; + + // open stream + std::ifstream os; + os.open(Rcpp::as (input_file),std::ifstream::in); + + // read file + SdpaFormatManager sdpaFormatManager; + sdpaFormatManager.loadSDPAFormatFile(os, _spectrahedron, c); + + std::vector matrices = _spectrahedron.getLMI().getMatrices(); + + // return spectrahedron and objective function + Rcpp::List _matrices; + + for (auto matrix : matrices) + _matrices.push_back(Rcpp::wrap(matrix)); + + Rcpp::List retList = Rcpp::List::create(Rcpp::Named("matrices") = _matrices , Rcpp::_["objFunction"] = Rcpp::wrap(c.getCoefficients())); + return Rcpp::wrap(retList); +} \ No newline at end of file diff --git a/src/ode_solve.cpp b/src/ode_solve.cpp index 07766d33..af20bed4 100644 --- a/src/ode_solve.cpp +++ b/src/ode_solve.cpp @@ -151,10 +151,10 @@ Rcpp::List ode_solve(Rcpp::Nullable n, Hpolytope HP(dim, Rcpp::as( Rcpp::as(Rcpp::as(domains) - [domain_name.c_str()]).field("A")), + [domain_name.c_str()]).slot("A")), Rcpp::as( Rcpp::as(Rcpp::as(domains) - [domain_name.c_str()]).field("b")) + [domain_name.c_str()]).slot("b")) ); HP.normalize(); diff --git a/src/rotating.cpp b/src/rotating.cpp index 607c0092..ebb9d115 100644 --- a/src/rotating.cpp +++ b/src/rotating.cpp @@ -30,7 +30,8 @@ //' //' @return A matrix that describes the rotated polytope // [[Rcpp::export]] -Rcpp::NumericMatrix rotating (Rcpp::Reference P, Rcpp::Nullable T = R_NilValue, +Rcpp::NumericMatrix rotating (Rcpp::Reference P, + Rcpp::Nullable T = R_NilValue, Rcpp::Nullable seed = R_NilValue){ typedef double NT; @@ -46,7 +47,25 @@ Rcpp::NumericMatrix rotating (Rcpp::Reference P, Rcpp::Nullable(P.slot("type")); + + if (type_str.compare(std::string("Hpolytope")) == 0) { + n = Rcpp::as(P.slot("A")).cols(); + type = 1; + } else if (type_str.compare(std::string("Vpolytope")) == 0) { + n = Rcpp::as(P.slot("V")).cols(); + type = 2; + } else if (type_str.compare(std::string("Zonotope")) == 0) { + n = Rcpp::as(P.slot("G")).cols(); + type = 3; + } else if (type_str.compare(std::string("VpolytopeIntersection")) == 0) { + throw Rcpp::exception("volesti does not support roatation of this kind of representation."); + } else { + throw Rcpp::exception("Unknown polytope representation!"); + } int seed_rcpp = (!seed.isNotNull()) ? std::chrono::system_clock::now() .time_since_epoch().count() @@ -55,7 +74,7 @@ Rcpp::NumericMatrix rotating (Rcpp::Reference P, Rcpp::Nullable(P.field("A")), Rcpp::as(P.field("b"))); + Hpolytope HP(n, Rcpp::as(P.slot("A")), Rcpp::as(P.slot("b"))); if (T.isNotNull()) { TransorfMat = Rcpp::as(T); HP.linear_transformIt(TransorfMat.inverse()); @@ -67,7 +86,7 @@ Rcpp::NumericMatrix rotating (Rcpp::Reference P, Rcpp::Nullable(P.field("V")), VT::Ones(Rcpp::as(P.field("V")).rows())); + Vpolytope VP(n, Rcpp::as(P.slot("V")), VT::Ones(Rcpp::as(P.slot("V")).rows())); if (T.isNotNull()) { TransorfMat = Rcpp::as(T); VP.linear_transformIt(TransorfMat.inverse()); @@ -79,7 +98,7 @@ Rcpp::NumericMatrix rotating (Rcpp::Reference P, Rcpp::Nullable(P.field("G")), VT::Ones(Rcpp::as(P.field("G")).rows())); + zonotope ZP(n, Rcpp::as(P.slot("G")), VT::Ones(Rcpp::as(P.slot("G")).rows())); if (T.isNotNull()) { TransorfMat = Rcpp::as(T); ZP.linear_transformIt(TransorfMat.inverse()); diff --git a/src/rounding.cpp b/src/rounding.cpp index 0b2cbe0f..5bb4977e 100644 --- a/src/rounding.cpp +++ b/src/rounding.cpp @@ -25,21 +25,21 @@ #include "preprocess/max_inscribed_ellipsoid_rounding.hpp" #include "extractMatPoly.h" -template +template < typename MT, - typename VT, - typename WalkType, - typename Polytope, - typename Point, - typename NT, + typename VT, + typename WalkType, + typename Polytope, + typename Point, + typename NT, typename RNGType > std::tuple apply_rounding(Polytope &P, std::string const& method_rcpp, unsigned int const& walkL, - std::pair &InnerBall, - RNGType &rng) + std::pair &InnerBall, + RNGType &rng) { std::tuple round_res; if (method_rcpp.compare(std::string("min_ellipsoid")) == 0) { @@ -62,7 +62,7 @@ std::tuple apply_rounding(Polytope &P, //' //' @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::List rounding (Rcpp::Reference P, Rcpp::Nullable method = R_NilValue, Rcpp::Nullable seed = R_NilValue){ @@ -76,7 +76,27 @@ Rcpp::List rounding (Rcpp::Reference P, typedef Eigen::Matrix VT; typedef Eigen::Matrix MT; - unsigned int n = P.field("dimension"), walkL = 2, type = P.field("type"); + unsigned int n; + unsigned int walkL=2; + unsigned int type; + + std::string type_str = Rcpp::as(P.slot("type")); + + if (type_str.compare(std::string("Hpolytope")) == 0) { + n = Rcpp::as(P.slot("A")).cols(); + type = 1; + } else if (type_str.compare(std::string("Vpolytope")) == 0) { + n = Rcpp::as(P.slot("V")).cols(); + type = 2; + } else if (type_str.compare(std::string("Zonotope")) == 0) { + n = Rcpp::as(P.slot("G")).cols(); + type = 3; + } else if (type_str.compare(std::string("VpolytopeIntersection")) == 0) { + throw Rcpp::exception("volesti does not support roatation of this kind of representation."); + } else { + throw Rcpp::exception("Unknown polytope representation!"); + } + std::string method_rcpp = std::string("isotropy"); if(method.isNotNull()) { method_rcpp = Rcpp::as(method); @@ -98,11 +118,10 @@ Rcpp::List rounding (Rcpp::Reference P, switch (type) { case 1: { // Hpolytope - - if (Rcpp::as(P.field("Aeq")).rows() > 0) { - throw Rcpp::exception("volesti supports rounding for full dimensional polytopes"); - } - Hpolytope HP(n, Rcpp::as(P.field("A")), Rcpp::as(P.field("b"))); + //if (Rcpp::as(P.slot("Aeq")).rows() > 0) { + // throw Rcpp::exception("volesti supports rounding for full dimensional polytopes"); + // } + Hpolytope HP(n, Rcpp::as(P.slot("A")), Rcpp::as(P.slot("b"))); InnerBall = HP.ComputeInnerBall(); if (InnerBall.second < 0.0) throw Rcpp::exception("Unable to compute a feasible point."); if (method_rcpp.compare(std::string("max_ellipsoid")) == 0) { @@ -115,7 +134,7 @@ Rcpp::List rounding (Rcpp::Reference P, } case 2: { // Vpolytope - Vpolytope VP(n, Rcpp::as(P.field("V")), VT::Ones(Rcpp::as(P.field("V")).rows())); + Vpolytope VP(n, Rcpp::as(P.slot("V")), VT::Ones(Rcpp::as(P.slot("V")).rows())); InnerBall = VP.ComputeInnerBall(); if (InnerBall.second < 0.0) throw Rcpp::exception("Unable to compute a feasible point."); round_res = apply_rounding(VP, method_rcpp, walkL, InnerBall, rng); @@ -124,7 +143,7 @@ Rcpp::List rounding (Rcpp::Reference P, } case 3: { // Zonotope - zonotope ZP(n, Rcpp::as(P.field("G")), VT::Ones(Rcpp::as(P.field("G")).rows())); + zonotope ZP(n, Rcpp::as(P.slot("G")), VT::Ones(Rcpp::as(P.slot("G")).rows())); InnerBall = ZP.ComputeInnerBall(); if (InnerBall.second < 0.0) throw Rcpp::exception("Unable to compute a feasible point."); round_res = apply_rounding(ZP, method_rcpp, walkL, InnerBall, rng); diff --git a/src/sample_points.cpp b/src/sample_points.cpp index d03c9b33..10f6dcee 100644 --- a/src/sample_points.cpp +++ b/src/sample_points.cpp @@ -2,7 +2,7 @@ // VolEsti (volume computation and sampling library) -// Copyright (c) 2012-2021 Vissarion Fisikopoulos +// Copyright (c) 2012-2024 Vissarion Fisikopoulos // Copyright (c) 2018-2021 Apostolos Chalkis //Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 and 2019 program. @@ -251,7 +251,7 @@ void sample_from_polytope(Polytope &P, int type, RNGType &rng, PointList &randPo //' \item{\code{BaW_rad} }{ The radius for the ball walk.} //' \item{\code{L} }{ The maximum length of the billiard trajectory or the radius for the step of dikin, vaidya or john walk.} //' \item{\code{solver} }{ Specify ODE solver for logconcave sampling. Options are i) leapfrog, ii) euler iii) runge-kutta iv) richardson} -//' \item{\code{step_size }{ Optionally chosen step size for logconcave sampling. Defaults to a theoretical value if not provided.} +//' \item{\code{step_size} }{ Optionally chosen step size for logconcave sampling. Defaults to a theoretical value if not provided.} //' } //' @param distribution Optional. A list that declares the target density and some related parameters as follows: //' \itemize{ @@ -293,7 +293,7 @@ void sample_from_polytope(Polytope &P, int type, RNGType &rng, PointList &randPo //' # gaussian distribution from the 2d unit simplex in H-representation with variance = 2 //' A = matrix(c(-1,0,0,-1,1,1), ncol=2, nrow=3, byrow=TRUE) //' b = c(0,0,1) -//' P = Hpolytope$new(A,b) +//' P = Hpolytope(A = A, b = b) //' points = sample_points(P, n = 100, distribution = list("density" = "gaussian", "variance" = 2)) //' //' # uniform points from the boundary of a 2-dimensional random H-polytope @@ -304,7 +304,7 @@ void sample_from_polytope(Polytope &P, int type, RNGType &rng, PointList &randPo //' //' @export // [[Rcpp::export]] -Rcpp::NumericMatrix sample_points(Rcpp::Nullable P, +Rcpp::NumericMatrix sample_points(Rcpp::Reference P, Rcpp::Nullable n, Rcpp::Nullable random_walk = R_NilValue, Rcpp::Nullable distribution = R_NilValue, @@ -321,15 +321,33 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P, typedef Eigen::Matrix VT; typedef Eigen::Matrix MT; - unsigned int type = Rcpp::as(P).field("type"), dim = Rcpp::as(P).field("dimension"), - walkL = 1, numpoints, nburns = 0; - RcppFunctor::GradientFunctor *F = NULL; RcppFunctor::FunctionFunctor *f = NULL; GaussianFunctor::GradientFunctor *G = NULL; GaussianFunctor::FunctionFunctor *g = NULL; bool functor_defined = true; + unsigned int dim; + unsigned int walkL = 1; + unsigned int type; + + std::string type_str = Rcpp::as(P.slot("type")); + + if (type_str.compare(std::string("Hpolytope")) == 0) { + dim = Rcpp::as(P.slot("A")).cols(); + type = 1; + } else if (type_str.compare(std::string("Vpolytope")) == 0) { + dim = Rcpp::as(P.slot("V")).cols(); + type = 2; + } else if (type_str.compare(std::string("Zonotope")) == 0) { + dim = Rcpp::as(P.slot("G")).cols(); + type = 3; + } else if (type_str.compare(std::string("VpolytopeIntersection")) == 0) { + dim = Rcpp::as(P.slot("V1")).cols(); + type = 4; + } else { + throw Rcpp::exception("Unknown polytope representation!"); + } RNGType rng(dim); if (seed.isNotNull()) { @@ -351,7 +369,7 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P, std::pair InnerBall; Point mode(dim); - numpoints = Rcpp::as(n); + unsigned int numpoints = Rcpp::as(n); if (numpoints <= 0) throw Rcpp::exception("The number of samples has to be a positive integer!"); if (!distribution.isNotNull() || !Rcpp::as(distribution).containsElementNamed("density")) { @@ -599,6 +617,7 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P, } } + unsigned int nburns = 0; if (Rcpp::as(random_walk).containsElementNamed("nburns")) { nburns = Rcpp::as(Rcpp::as(random_walk)["nburns"]); if (nburns < 0) { @@ -609,8 +628,8 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P, switch(type) { case 1: { // Hpolytope - Hpolytope HP(dim, Rcpp::as(Rcpp::as(P).field("A")), - Rcpp::as(Rcpp::as(P).field("b"))); + Hpolytope HP(dim, Rcpp::as(P.slot("A")), + Rcpp::as(P.slot("b"))); InnerBall = HP.ComputeInnerBall(); if (InnerBall.second < 0.0) throw Rcpp::exception("Unable to compute a feasible point."); @@ -637,8 +656,8 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P, } case 2: { // Vpolytope - Vpolytope VP(dim, Rcpp::as(Rcpp::as(P).field("V")), - VT::Ones(Rcpp::as(Rcpp::as(P).field("V")).rows())); + Vpolytope VP(dim, Rcpp::as(P.slot("V")), + VT::Ones(Rcpp::as(P.slot("V")).rows())); InnerBall = VP.ComputeInnerBall(); if (InnerBall.second < 0.0) throw Rcpp::exception("Unable to compute a feasible point."); @@ -658,8 +677,8 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P, } case 3: { // Zonotope - zonotope ZP(dim, Rcpp::as(Rcpp::as(P).field("G")), - VT::Ones(Rcpp::as(Rcpp::as(P).field("G")).rows())); + zonotope ZP(dim, Rcpp::as(P.slot("G")), + VT::Ones(Rcpp::as(P.slot("G")).rows())); InnerBall = ZP.ComputeInnerBall(); if (InnerBall.second < 0.0) throw Rcpp::exception("Unable to compute a feasible point."); @@ -679,10 +698,10 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P, } case 4: { // Intersection of two V-polytopes - Vpolytope VP1(dim, Rcpp::as(Rcpp::as(P).field("V1")), - VT::Ones(Rcpp::as(Rcpp::as(P).field("V1")).rows())); - Vpolytope VP2(dim, Rcpp::as(Rcpp::as(P).field("V2")), - VT::Ones(Rcpp::as(Rcpp::as(P).field("V2")).rows())); + Vpolytope VP1(dim, Rcpp::as(P.slot("V1")), + VT::Ones(Rcpp::as(P.slot("V1")).rows())); + Vpolytope VP2(dim, Rcpp::as(P.slot("V2")), + VT::Ones(Rcpp::as(P.slot("V2")).rows())); InterVP VPcVP(VP1, VP2); if (!VPcVP.is_feasible()) throw Rcpp::exception("Empty set!"); diff --git a/src/spectrahedron.cpp b/src/spectrahedron.cpp index e7926597..f084fd65 100644 --- a/src/spectrahedron.cpp +++ b/src/spectrahedron.cpp @@ -35,13 +35,13 @@ //' 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$new(lmi); +//' S = Spectrahedron(matrices = lmi) //' objFunction = c(1,1) //' writeSdpaFormatFile(S, objFunction, "output.txt") //' } //' @export // [[Rcpp::export]] -void writeSdpaFormatFile(Rcpp::Nullable spectrahedron = R_NilValue, +void writeSdpaFormatFile(Rcpp::Reference spectrahedron = R_NilValue, Rcpp::Nullable objectiveFunction = R_NilValue, Rcpp::Nullable outputFile = R_NilValue) { @@ -54,7 +54,7 @@ void writeSdpaFormatFile(Rcpp::Nullable spectrahedron = R_NilVa typedef LMI LMI; typedef Spectrahedron Spectrahedron; - std::vector matrices = Rcpp::as > (Rcpp::as (spectrahedron).field("matrices")); + std::vector matrices = Rcpp::as > (Rcpp::as (spectrahedron).slot("matrices")); LMI lmi(matrices); Spectrahedron _spectrahedron(lmi); Point c(Rcpp::as (objectiveFunction)); diff --git a/src/volume.cpp b/src/volume.cpp index 0a47210a..49dd9960 100644 --- a/src/volume.cpp +++ b/src/volume.cpp @@ -2,7 +2,7 @@ // VolEsti (volume computation and sampling library) -// Copyright (c) 2012-2020 Vissarion Fisikopoulos +// Copyright (c) 2012-2024 Vissarion Fisikopoulos // Copyright (c) 2018-2020 Apostolos Chalkis //Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 and 2019 program. @@ -29,8 +29,8 @@ enum rounding_type {none, min_ellipsoid, max_ellipsoid, isotropy}; template std::pair generic_volume(Polytope& P, RNGType &rng, unsigned int walk_length, NT e, - volume_algorithms const& algo, unsigned int win_len, - rounding_type const& rounding, random_walks const& walk) + volume_algorithms const& algo, unsigned int win_len, + rounding_type const& rounding, random_walks const& walk) { typedef typename Polytope::MT MT; typedef typename Polytope::VT VT; @@ -181,6 +181,7 @@ std::pair generic_volume(Polytope& P, RNGType &rng, unsigned int //' \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{250} for CB with BiW and \eqn{400+3d^2} for CB and any other random walk and \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 Optional. A string parameter to request a rounding method to be applied in the input polytope before volume computation: a) \code{'min_ellipsoid'}, b) \code{'svd'}, c) \code{'max_ellipsoid'} and d) \code{'none'} for no rounding. //' @param seed Optional. A fixed seed for the number generator. @@ -210,9 +211,8 @@ std::pair generic_volume(Polytope& P, RNGType &rng, unsigned int //' @export // [[Rcpp::export]] Rcpp::List volume (Rcpp::Reference P, - Rcpp::Nullable settings = R_NilValue, - Rcpp::Nullable rounding = R_NilValue, - Rcpp::Nullable seed = R_NilValue) { + Rcpp::Nullable settings = R_NilValue, + Rcpp::Nullable rounding = R_NilValue) { typedef double NT; typedef Cartesian Kernel; @@ -224,12 +224,30 @@ Rcpp::List volume (Rcpp::Reference P, typedef IntersectionOfVpoly InterVP; typedef Eigen::Matrix VT; typedef Eigen::Matrix MT; - unsigned int n = P.field("dimension"), walkL, type = P.field("type"); + + unsigned int n, walkL, type; + std::string type_str = Rcpp::as(P.slot("type")); + + if (type_str.compare(std::string("Hpolytope")) == 0) { + n = Rcpp::as(P.slot("A")).cols(); + type = 1; + } else if (type_str.compare(std::string("Vpolytope")) == 0) { + n = Rcpp::as(P.slot("V")).cols(); + type = 2; + } else if (type_str.compare(std::string("Zonotope")) == 0) { + n = Rcpp::as(P.slot("G")).cols(); + type = 3; + } else if (type_str.compare(std::string("VpolytopeIntersection")) == 0) { + n = Rcpp::as(P.slot("V1")).cols(); + type = 4; + } else { + throw Rcpp::exception("Unknown polytope representation!"); + } RNGType rng(n); - if (seed.isNotNull()) { - unsigned seed_rcpp = Rcpp::as(seed); - rng.set_seed(seed_rcpp); + if (Rcpp::as(settings).containsElementNamed("seed")) { + unsigned seed_tmp = Rcpp::as(Rcpp::as(settings)["seed"]); + rng.set_seed(seed_tmp); } bool round = false, hpoly = false; @@ -333,19 +351,19 @@ Rcpp::List volume (Rcpp::Reference P, switch(type) { case 1: { // Hpolytope - Hpolytope HP(n, Rcpp::as(P.field("A")), Rcpp::as(P.field("b"))); + Hpolytope HP(n, Rcpp::as(P.slot("A")), Rcpp::as(P.slot("b"))); pair_vol = generic_volume(HP, rng, walkL, e, algo, win_len, rounding_method, walk); break; } case 2: { // Vpolytope - Vpolytope VP(n, Rcpp::as(P.field("V")), VT::Ones(Rcpp::as(P.field("V")).rows())); + Vpolytope VP(n, Rcpp::as(P.slot("V")), VT::Ones(Rcpp::as(P.slot("V")).rows())); pair_vol = generic_volume(VP, rng, walkL, e, algo, win_len, rounding_method, walk); break; } case 3: { // Zonotope - zonotope ZP(n, Rcpp::as(P.field("G")), VT::Ones(Rcpp::as(P.field("G")).rows())); + zonotope ZP(n, Rcpp::as(P.slot("G")), VT::Ones(Rcpp::as(P.slot("G")).rows())); if (Rcpp::as(settings).containsElementNamed("hpoly")) { hpoly = Rcpp::as(Rcpp::as(settings)["hpoly"]); if (hpoly && (algo == CG || algo == SOB)) @@ -388,14 +406,15 @@ Rcpp::List volume (Rcpp::Reference P, } case 4: { // Intersection of two V-polytopes - Vpolytope VP1(n, Rcpp::as(P.field("V1")), VT::Ones(Rcpp::as(P.field("V1")).rows())); - Vpolytope VP2(n, Rcpp::as(P.field("V2")), VT::Ones(Rcpp::as(P.field("V2")).rows())); + Vpolytope VP1(n, Rcpp::as(P.slot("V1")), VT::Ones(Rcpp::as(P.slot("V1")).rows())); + Vpolytope VP2(n, Rcpp::as(P.slot("V2")), VT::Ones(Rcpp::as(P.slot("V2")).rows())); InterVP VPcVP; - if (!seed.isNotNull()) { - InterVP VPcVP(VP1, VP2); + if (Rcpp::as(settings).containsElementNamed("seed")) { + unsigned seed_tmp = Rcpp::as(Rcpp::as(settings)["seed"]); + rng.set_seed(seed_tmp); + InterVP VPcVP(VP1, VP2, seed_tmp); } else { - unsigned seed3 = Rcpp::as(seed); - InterVP VPcVP(VP1, VP2, seed3); + InterVP VPcVP(VP1, VP2); } if (!VPcVP.is_feasible()) throw Rcpp::exception("Empty set!"); pair_vol = generic_volume(VPcVP, rng, walkL, e, algo, win_len, rounding_method, walk); diff --git a/src/write_sdpa_format_file.cpp b/src/write_sdpa_format_file.cpp new file mode 100644 index 00000000..0be3f969 --- /dev/null +++ b/src/write_sdpa_format_file.cpp @@ -0,0 +1,69 @@ +// VolEsti (volume computation and sampling library) + +// Copyright (c) 2012-2020 Vissarion Fisikopoulos +// Copyright (c) 2020 Apostolos Chalkis + +//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. + +// Licensed under GNU LGPL.3, see LICENCE file + +#include +#include +#include +#include "cartesian_geom/cartesian_kernel.h" +#include +#include +#include +#include +#include "convex_bodies/spectrahedra/LMI.h" +#include "convex_bodies/spectrahedra/spectrahedron.h" +#include "SDPAFormatManager.h" + +//' 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 +// [[Rcpp::export]] +void write_sdpa_format_file(Rcpp::Reference spectrahedron, + Rcpp::NumericVector objective_function, + std::string output_file) { + + typedef double NT; + typedef Eigen::Matrix VT; + typedef Eigen::Matrix MT; + typedef Cartesian Kernel; + typedef typename Kernel::Point Point; + typedef boost::mt19937 RNGType; + typedef LMI LMI; + typedef Spectrahedron Spectrahedron; + + std::vector matrices = Rcpp::as > (spectrahedron.slot("matrices")); + LMI lmi(matrices); + Spectrahedron _spectrahedron(lmi); + Point c(Rcpp::as (objective_function)); + + std::filebuf fb; + fb.open(output_file, std::ios::out); + std::ostream os(&fb); + + SdpaFormatManager sdpaFormatManager; + sdpaFormatManager.writeSDPAFormatFile(os, _spectrahedron, c); + + return; +} \ No newline at end of file diff --git a/src/zonotope_approximation.cpp b/src/zonotope_approximation.cpp index 237ec96e..7debd237 100644 --- a/src/zonotope_approximation.cpp +++ b/src/zonotope_approximation.cpp @@ -42,7 +42,17 @@ Rcpp::List zono_approx (Rcpp::Reference Z, typedef Zonotope zonotope; typedef Eigen::Matrix VT; typedef Eigen::Matrix MT; - int n = Rcpp::as(Z.field("dimension")), k = Rcpp::as(Z.field("G")).rows(), win_len = 200, walkL = 1; + + int k = Rcpp::as(Z.slot("G")).rows(), win_len = 250, walkL = 1; + + std::string type = Rcpp::as(Z.slot("type")); + + int n; + if (type.compare(std::string("Zonotope")) == 0) { + n = Rcpp::as(Z.slot("G")).cols(); + } else { + throw Rcpp::exception("This is not a zonotope."); + } RNGType rng(n); if (seed.isNotNull()) { @@ -54,10 +64,10 @@ Rcpp::List zono_approx (Rcpp::Reference Z, bool hpoly = false; MT X(n, 2 * k); - X << Rcpp::as(Z.field("G")).transpose(), -Rcpp::as(Z.field("G")).transpose(); + X << Rcpp::as(Z.slot("G")).transpose(), -Rcpp::as(Z.slot("G")).transpose(); Eigen::JacobiSVD svd(X * X.transpose(), Eigen::ComputeFullU | Eigen::ComputeFullV); MT G(k, 2 * n); - G << Rcpp::as(Z.field("G")) * svd.matrixU(), Rcpp::as(Z.field("G")) * svd.matrixU(); + G << Rcpp::as(Z.slot("G")) * svd.matrixU(), Rcpp::as(Z.slot("G")) * svd.matrixU(); VT Gred_ii = G.transpose().cwiseAbs().rowwise().sum(); MT A(n, 2 * n); A << -MT::Identity(n, n), MT::Identity(n, n); @@ -79,7 +89,7 @@ Rcpp::List zono_approx (Rcpp::Reference Z, win_len = (!Rcpp::as(settings).containsElementNamed("win_len")) ? 200 : Rcpp::as( Rcpp::as(settings)["win_len"]); - zonotope ZP(n, Rcpp::as(Z.field("G")), VT::Ones(Rcpp::as(Z.field("G")).rows())); + zonotope ZP(n, Rcpp::as(Z.slot("G")), VT::Ones(Rcpp::as(Z.slot("G")).rows())); if (Rcpp::as(settings).containsElementNamed("hpoly")) { hpoly = Rcpp::as(Rcpp::as(settings)["hpoly"]);