Make R work with R classes and changes from volesti develop

.github/workflows/R-CMD-check-macOS.yml

@@ -37,16 +37,13 @@ jobs:
 
       - uses: r-lib/actions/setup-pandoc@v2
 
-      - name: Fetch and update submodule
-        run: git submodule update --init --recursive
-
       - name: Install dependencies
         run: Rscript -e "install.packages(c('devtools', dependencies=TRUE))" -e "install.packages(c('rcmdcheck', 'devtools', 'Rcpp', 'RcppEigen', 'BH', 'testthat', 'downloader', 'xfun'))";
 
       - name: Check
         env:
           _R_CHECK_CRAN_INCOMING_REMOTE_: false
-        run: Rscript -e "library(rcmdcheck)" -e "rcmdcheck::rcmdcheck(args = c('--no-manual', '--no-install'), error_on = 'error', check_dir = 'check')"
+        run: Rscript -e "library(rcmdcheck)" -e "rcmdcheck::rcmdcheck(args = c('--no-manual'), error_on = 'warning', check_dir = 'check')"
 
       - name: Upload check results
         if: failure()

.github/workflows/R-CMD-check-ubuntu.yml

@@ -39,16 +39,13 @@ jobs:
 
       - uses: r-lib/actions/setup-pandoc@v2
 
-      - name: Fetch and update submodule
-        run: git submodule update --init --recursive
-
       - name: Install dependencies
         run: Rscript -e "install.packages(c('testthat', 'pkgload', 'rcmdcheck', 'devtools', 'Rcpp', 'RcppEigen', 'BH', 'downloader', 'xfun', dependencies=TRUE))";
 
       - name: Check
         env:
           _R_CHECK_CRAN_INCOMING_REMOTE_: false
-        run: Rscript -e "library(rcmdcheck)" -e "rcmdcheck::rcmdcheck(args = c('--no-manual', '--no-install'), error_on = 'error', check_dir = 'check')"
+        run: Rscript -e "library(rcmdcheck)" -e "rcmdcheck::rcmdcheck(args = c('--no-manual'), error_on = 'warning', check_dir = 'check')"
 
       - name: Upload check results
         if: failure()

.github/workflows/R-CMD-check-windows.yml

@@ -36,16 +36,13 @@ jobs:
 
       - uses: r-lib/actions/setup-pandoc@v2
 
-      - name: Fetch and update submodule
-        run: git submodule update --init --recursive
-
       - name: Install dependencies
         run: Rscript -e "install.packages(c('devtools', dependencies=TRUE))" -e "install.packages(c('rcmdcheck', 'devtools', 'Rcpp', 'RcppEigen', 'BH', 'testthat', 'downloader', 'xfun'))"
 
       - name: Check
         env:
           _R_CHECK_CRAN_INCOMING_REMOTE_: false
-        run: Rscript -e "library(rcmdcheck)" -e "rcmdcheck::rcmdcheck(args = c('--no-manual', '--no-install'), error_on = 'error', check_dir = 'check')"
+        run: Rscript -e "library(rcmdcheck)" -e "rcmdcheck::rcmdcheck(args = c('--no-manual'), error_on = 'warning', check_dir = 'check')"
 
       - name: Upload check results
         if: failure()

DESCRIPTION

@@ -12,13 +12,13 @@ Description: Provides an R interface for 'volesti' C++ package. 'volesti' comput
              for sampling, rounding and rotating polytopes. Moreover, 'volesti' provides algorithms for
              estimating copulas useful in computational finance. Methods implemented in 'volesti' are described
              in A. Chalkis and V. Fisikopoulos (2022) <doi:10.32614/RJ-2021-077> and references therein.
-Version: 1.1.2-6
-Date: 2023-04-11
+Version: 1.2.0
+Date: 2024-02-29
 Maintainer: Vissarion Fisikopoulos <[email protected]>
 Depends: Rcpp (>= 0.12.17)
 Imports: methods, stats
 LinkingTo: Rcpp, RcppEigen, BH
 Suggests: testthat
 Encoding: UTF-8
-RoxygenNote: 7.2.3
-BugReports: https://github.com/GeomScale/volesti/issues
+RoxygenNote: 7.3.1
+BugReports: https://github.com/GeomScale/Rvolesti/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)

NEWS.md

@@ -50,3 +50,15 @@
 # volesti 1.1.2-6
 
 - Fix UBSAN issues (lp_presolve)
+
+# volesti 1.2.0
+
+- New functions: dinvweibull_with_loc, ess, estimtate_lipschitz_constant, gen_birkhoff, geweke
+ode_solve, pinvweibull_with_loc, psrf_multivariate, psrf_univariate, raftery
+
+- New features in sample_points function:
+ a) new walks: i) Dikin walk, ii) Vaidya walk, iii) John walk,
+  iv) Hamiltonian Monte Carlo (HMC) for general logconcave densities
+  v) Underdamped Langevin Dynamics (ULD) for general logconcave densities (using the Randomized Midpoint Method)
+  vi) Exact Hamiltonian Monte Carlo with reflections (spherical Gaussian or exponential distribution)
+ b) new distributions: i) exponential, ii) general logconcave

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.
@@ -191,7 +202,9 @@ inner_ball <- function(P, lpsolve = NULL) {
 #' @return A list which contains elements "x_1", ..., "x_n" representing each derivative results. Each "x_i" corresponds to a d x n matrix where each column represents a certain timestep of the solver.
 #'
 #' @examples
-#' # Please visit the examples directory on examples demonstrating usage of the ODE solvers.
+#' F <- function (x) (-x)
+#' initial_conditions <- list("x_1" = c(0), "x_2" = c(1))
+#' states <- ode_solve(dimension=1, n=1000, F=F, initial_time=0, step_size=0.01, order=2, method="leapfrog", initial_conditions=initial_conditions, domains = list())
 #'
 #' @export
 ode_solve <- function(n, step_size, order, dimension, initial_time, F, method, domains = NULL, initial_conditions = NULL) {
@@ -305,7 +318,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 +360,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 +389,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,9 +425,9 @@ 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.
 #'
 #' @references \cite{I.Z.Emiris and V. Fisikopoulos,
 #' \dQuote{Practical polytope volume approximation,} \emph{ACM Trans. Math. Soft.,} 2018.},
@@ -439,8 +452,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))
-  
-}
+
+}