Make R work with R classes and changes from develop

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

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)

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

R/compute_indicators.R

@@ -1,58 +1,87 @@
 #' 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) {
         compRet[j] = compRet[j] * (1 + returns[k, j])
       }
       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))
-  
-}
+
+}