From e859479fdafb4bb733fdae97c49c6489818cd200 Mon Sep 17 00:00:00 2001 From: Tolis Date: Wed, 12 Feb 2020 20:33:19 +0200 Subject: [PATCH 01/13] merge two copulas functions into one. --- R-proj/NAMESPACE | 3 +- R-proj/src/{copula1.cpp => copula.cpp} | 52 ++++++++++++++++++++++---- 2 files changed, 45 insertions(+), 10 deletions(-) rename R-proj/src/{copula1.cpp => copula.cpp} (50%) diff --git a/R-proj/NAMESPACE b/R-proj/NAMESPACE index 7fa456964..8699d7b67 100644 --- a/R-proj/NAMESPACE +++ b/R-proj/NAMESPACE @@ -1,7 +1,6 @@ # Generated by roxygen2: do not edit by hand -export(copula1) -export(copula2) +export(copula) export(exact_vol) export(file_to_polytope) export(frustum_of_simplex) diff --git a/R-proj/src/copula1.cpp b/R-proj/src/copula.cpp similarity index 50% rename from R-proj/src/copula1.cpp rename to R-proj/src/copula.cpp index 0093d81b8..053f69660 100644 --- a/R-proj/src/copula1.cpp +++ b/R-proj/src/copula.cpp @@ -9,6 +9,7 @@ #include #include +#include "ellipsoids.h" #include "cartesian_geom/cartesian_kernel.h" #include #include @@ -17,10 +18,11 @@ //' Construct a copula using uniform sampling from the unit simplex //' -//' Given two families of parallel hyperplanes intersecting the canonical simplex, this function uniformly samples from the canonical simplex and construct an approximation of the bivariate probability distribution, called copula. +//' Given two families of parallel hyperplanes or a family of parallel hyperplanes and a family of concentric ellispoids centered at the origin intersecting the canonical simplex, this function uniformly samples from the canonical simplex and construct an approximation of the bivariate probability distribution, called copula. //' //' @param h1 A \eqn{d}-dimensional vector that describes the direction of the first family of parallel hyperplanes. -//' @param h2 A \eqn{d}-dimensional vector that describes the direction of the second family of parallel hyperplanes. +//' @param h2 Optional. A \eqn{d}-dimensional vector that describes the direction of the second family of parallel hyperplanes. +//' @param E Optional. The \eqn{d\times d} symmetric positive semidefine matrix that describes the family of concentric ellipsoids centered at the origin. //' @param numSlices The number of the slices for the copula. Default value is 100. //' @param N The number of points to sample. Default value is \eqn{4\cdot 10^6}. //' @@ -34,16 +36,26 @@ //' h1 = h1 / 1000 //' h2=runif(n = 10, min = 1, max = 1000) //' h2 = h2 / 1000 -//' cop = copula1(h1=h1, h2=h2, numSlices = 10, N = 100000) +//' cop = copula(h1=h1, h2=h2, numSlices = 10, N = 100000) +//' +//' # compute a copula for a family of parallel hyperplanes and a family of conentric ellipsoids +//' h = runif(n = 10, min = 1, max = 1000) +//' h = h / 1000 +//' E = replicate(10, rnorm(20)) +//' E = cov(E) +//' cop = copula(h1=h, E=E, numSlices=10, N=100000) +//' //' @export // [[Rcpp::export]] -Rcpp::NumericMatrix copula1 (Rcpp::NumericVector h1, Rcpp::NumericVector h2, - Rcpp::Nullable numSlices, Rcpp::Nullable N){ +Rcpp::NumericMatrix copula (Rcpp::NumericVector h1 = R_NilValue, Rcpp::NumericVector h2 = R_NilValue, + Rcpp::NumericMatrix E = R_NilValue, Rcpp::Nullable numSlices = R_NilValue, + Rcpp::Nullable N = R_NilValue){ typedef double NT; typedef Cartesian Kernel; typedef typename Kernel::Point Point; typedef boost::mt19937 RNGType; + typedef Eigen::Matrix MT; unsigned int num_slices = 100, numpoints = 4000000; if (numSlices.isNotNull()) { @@ -56,12 +68,36 @@ Rcpp::NumericMatrix copula1 (Rcpp::NumericVector h1, Rcpp::NumericVector h2, Rcpp::NumericMatrix copula(num_slices, num_slices); std::vector > StdCopula; - unsigned int dim = h1.size(), i, j; + unsigned int dim = Rcpp::as >(h1).size(), i, j; + + if(!h1.isNotNull()) { + + throw Rcpp::exception("You have to give at least one normal of a hyperplane!"); + + } std::vector hyp1 = Rcpp::as >(h1); - std::vector hyp2 = Rcpp::as >(h2); + if (h2.isNotNull()) { - StdCopula = twoParHypFam(dim, numpoints, num_slices, hyp1, hyp2); + std::vector hyp2 = Rcpp::as < std::vector < NT > > (h2); + StdCopula = twoParHypFam(dim, numpoints, num_slices, hyp1, hyp2); + + } else if (E.isNotNull()) { + + std::vector > Gin(dim, std::vector(dim)); + MT EE = Rcpp::as(E); + for (i=0; i(dim, numpoints, num_slices, hyp1, Ell); + } else { + + throw Rcpp::exception("Wrong inputs"); + + } for(i=0; i Date: Wed, 12 Feb 2020 20:47:07 +0200 Subject: [PATCH 02/13] add compute_indicators R function. --- R-proj/NAMESPACE | 1 + R-proj/R/compute_indicators.R | 73 +++++++++++++++++++++++++++++++++++ 2 files changed, 74 insertions(+) create mode 100644 R-proj/R/compute_indicators.R diff --git a/R-proj/NAMESPACE b/R-proj/NAMESPACE index 8699d7b67..9710f5374 100644 --- a/R-proj/NAMESPACE +++ b/R-proj/NAMESPACE @@ -1,6 +1,7 @@ # Generated by roxygen2: do not edit by hand export(copula) +export(compute_indicators) export(exact_vol) export(file_to_polytope) export(frustum_of_simplex) diff --git a/R-proj/R/compute_indicators.R b/R-proj/R/compute_indicators.R new file mode 100644 index 000000000..37fdfbdec --- /dev/null +++ b/R-proj/R/compute_indicators.R @@ -0,0 +1,73 @@ +compute_indicators <- function(MatReturns, win_len = NULL, numSlices = NULL, n = NULL) { + + if (is.null(win_len)) win_len = 60 + if (is.null(numSlices)) numSlices = 100 + if (is.null(N)) N = 500000 + + nrows = dim(MatReturns)[1] + nassets = dim(MatReturns)[2] + wl = win_len-1 + + indicators = c() + print(nrows-wl) + for (i in 1:(nrows-wl)) { + + W=i:(i+wl) + #nR = MatReturns[W,] + E = cov(MatReturns[W,]) + + compRet = rep(1,nassets) + for (j in 1:nassets) { + for (k in W) { + compRet[j] = compRet[j] * (1 + MatReturns[k, j]) + } + compRet[j] = compRet[j] - 1 + } + + cop = copula2(h = compRet, E = E, numSlices = numSlices, N = N) + blue_mass = 0 + red_mass = 0 + + for (row in 1:100) { + for (col in 1:100) { + if (row-col<=0.2*numSlices && row-col>=-0.2*numSlices) { + if (row+col<0.8*numSlices || row+col>1.2*numSlices) { + red_mass = red_mass + cop[row,col] + } + } else { + if (row+col>=0.8*numSlices && row+col<=1.2*numSlices) { + blue_mass = blue_mass + cop[row,col] + } + } + } + } + indicators = c(indicators, blue_mass / red_mass) + print(length(indicators)) + } + + n = length(indicators) + + index = 0 + set_index = FALSE + col = rep("normal",n) + + for (i in 1:n) { + + if(indicators[i]>1 && !set_index){ + index = i + set_index = TRUE + } else if (indicators[i]<1) { + if(set_index){ + if(i-index+1>30 && i-index+1<60){ + col[index:i] = "warning" + } else if(i-index+1>60) { + col[index:i] = "crisis" + } + } + set_index = FALSE + } + } + + return(list("indicators" = indicatos, market_states = col) + +} From 753975228d936af7416c8e761a7fb093b7fc9a7c Mon Sep 17 00:00:00 2001 From: Tolis Date: Wed, 12 Feb 2020 21:19:11 +0200 Subject: [PATCH 03/13] develop compute_indicator R function and copula.cpp in R-proj --- R-proj/DESCRIPTION | 2 +- R-proj/R/compute_indicators.R | 46 ++++++++++++++++++++++++----------- R-proj/src/copula.cpp | 6 ++--- 3 files changed, 36 insertions(+), 18 deletions(-) diff --git a/R-proj/DESCRIPTION b/R-proj/DESCRIPTION index ab21adf38..d3188e83a 100644 --- a/R-proj/DESCRIPTION +++ b/R-proj/DESCRIPTION @@ -17,7 +17,7 @@ Maintainer: Vissarion Fisikopoulos Depends: Rcpp (>= 0.12.17) Imports: methods LinkingTo: Rcpp, RcppEigen, BH -Suggests: testthat +Suggests: testthat, tawny Encoding: UTF-8 RoxygenNote: 6.1.1 BugReports: https://github.com/GeomScale/volume_approximation/issues diff --git a/R-proj/R/compute_indicators.R b/R-proj/R/compute_indicators.R index 37fdfbdec..817a8f069 100644 --- a/R-proj/R/compute_indicators.R +++ b/R-proj/R/compute_indicators.R @@ -1,41 +1,59 @@ -compute_indicators <- function(MatReturns, win_len = NULL, numSlices = NULL, n = NULL) { +#' Compute an indicator for each time period that describes the state of a market. +#' +#' Given a matrix that contains row-wise the assets' returns and a sliding window W, this function computes an approximation of the joint distribution (copula) between portfolios' return and volatility in each time period implied by W. For each copula it computes an indicator: large value corresponds to a crisis period and a small value to a normal period. If 60 consecutive indicators are larger than 1. +#' +#' @param MatReturns A \eqn{d}-dimensional vector that describes the direction of the first family of parallel hyperplanes. +#' @param W 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}. +#' +#' @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. +#' +#' @export +compute_indicators <- function(MatReturns, W = NULL, M = NULL, N = NULL) { - if (is.null(win_len)) win_len = 60 - if (is.null(numSlices)) numSlices = 100 + if (is.null(W)) W = 60 + if (is.null(M)) M = 100 if (is.null(N)) N = 500000 nrows = dim(MatReturns)[1] nassets = dim(MatReturns)[2] - wl = win_len-1 + wl = W-1 indicators = c() print(nrows-wl) for (i in 1:(nrows-wl)) { - W=i:(i+wl) - #nR = MatReturns[W,] - E = cov(MatReturns[W,]) + Win=i:(i+wl) + if("tawny" %in% rownames(installed.packages()) == FALSE) { + E = cov(MatReturns[Win,]) + } else { + E = tawny::cov.shrink(MatReturns[Win,]) + } compRet = rep(1,nassets) for (j in 1:nassets) { - for (k in W) { + for (k in Win) { compRet[j] = compRet[j] * (1 + MatReturns[k, j]) } compRet[j] = compRet[j] - 1 } - cop = copula2(h = compRet, E = E, numSlices = numSlices, N = N) + cop = copula(h1 = compRet, E = E, numSlices = M, N = N) blue_mass = 0 red_mass = 0 - for (row in 1:100) { - for (col in 1:100) { - if (row-col<=0.2*numSlices && row-col>=-0.2*numSlices) { - if (row+col<0.8*numSlices || row+col>1.2*numSlices) { + for (row in 1:M) { + for (col in 1:M) { + if (row-col<=0.2*M && row-col>=-0.2*M) { + if (row+col<0.8*M || row+col>1.2*M) { red_mass = red_mass + cop[row,col] } } else { - if (row+col>=0.8*numSlices && row+col<=1.2*numSlices) { + if (row+col>=0.8*M && row+col<=1.2*M) { blue_mass = blue_mass + cop[row,col] } } diff --git a/R-proj/src/copula.cpp b/R-proj/src/copula.cpp index 053f69660..80d271bd3 100644 --- a/R-proj/src/copula.cpp +++ b/R-proj/src/copula.cpp @@ -23,8 +23,8 @@ //' @param h1 A \eqn{d}-dimensional vector that describes the direction of the first family of parallel hyperplanes. //' @param h2 Optional. A \eqn{d}-dimensional vector that describes the direction of the second family of parallel hyperplanes. //' @param E Optional. The \eqn{d\times d} symmetric positive semidefine matrix that describes the family of concentric ellipsoids centered at the origin. -//' @param numSlices The number of the slices for the copula. Default value is 100. -//' @param N The number of points to sample. Default value is \eqn{4\cdot 10^6}. +//' @param numSlices The number of the slices for the copula. The default value is 100. +//' @param N The number of points to sample. The default value is \eqn{5\cdot 10^5}. //' //' @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.} @@ -56,7 +56,7 @@ Rcpp::NumericMatrix copula (Rcpp::NumericVector h1 = R_NilValue, Rcpp::NumericVe typedef typename Kernel::Point Point; typedef boost::mt19937 RNGType; typedef Eigen::Matrix MT; - unsigned int num_slices = 100, numpoints = 4000000; + unsigned int num_slices = 100, numpoints = 500000; if (numSlices.isNotNull()) { num_slices = Rcpp::as(numSlices); From e6cd83f5e168aa455d7a41e2ed45a2b3a90a9b45 Mon Sep 17 00:00:00 2001 From: Tolis Date: Wed, 12 Feb 2020 21:28:03 +0200 Subject: [PATCH 04/13] improve roxygen comments in compute_indicator R function and copula.cpp in R-proj --- R-proj/R/compute_indicators.R | 6 ++++-- R-proj/src/copula.cpp | 26 +++++++++++++------------- 2 files changed, 17 insertions(+), 15 deletions(-) diff --git a/R-proj/R/compute_indicators.R b/R-proj/R/compute_indicators.R index 817a8f069..eaac88c69 100644 --- a/R-proj/R/compute_indicators.R +++ b/R-proj/R/compute_indicators.R @@ -1,6 +1,8 @@ #' Compute an indicator for each time period that describes the state of a market. #' -#' Given a matrix that contains row-wise the assets' returns and a sliding window W, this function computes an approximation of the joint distribution (copula) between portfolios' return and volatility in each time period implied by W. For each copula it computes an indicator: large value corresponds to a crisis period and a small value to a normal period. If 60 consecutive indicators are larger than 1. +#' Given a matrix that contains row-wise the assets' returns and a sliding window W, this function computes an approximation of the joint distribution (copula) between portfolios' return and volatility in each time period implied by W. +#' For each copula it computes an indicator: large value corresponds to a crisis period and a small value to a normal period. +#' The periods over which the indicator is greater than 1 for more than 60 consecutives sliding windows are warnings and for more than 100 are crisis. The sliding window is shifted by one day. #' #' @param MatReturns A \eqn{d}-dimensional vector that describes the direction of the first family of parallel hyperplanes. #' @param W Optional. The length of the sliding window. The default value is 60. @@ -42,7 +44,7 @@ compute_indicators <- function(MatReturns, W = NULL, M = NULL, N = NULL) { compRet[j] = compRet[j] - 1 } - cop = copula(h1 = compRet, E = E, numSlices = M, N = N) + cop = copula(R1 = compRet, Sigma = E, M = M, N = N) blue_mass = 0 red_mass = 0 diff --git a/R-proj/src/copula.cpp b/R-proj/src/copula.cpp index 80d271bd3..b79866d69 100644 --- a/R-proj/src/copula.cpp +++ b/R-proj/src/copula.cpp @@ -23,7 +23,7 @@ //' @param h1 A \eqn{d}-dimensional vector that describes the direction of the first family of parallel hyperplanes. //' @param h2 Optional. A \eqn{d}-dimensional vector that describes the direction of the second family of parallel hyperplanes. //' @param E Optional. The \eqn{d\times d} symmetric positive semidefine matrix that describes the family of concentric ellipsoids centered at the origin. -//' @param numSlices The number of the slices for the copula. The default value is 100. +//' @param M The number of the slices for the copula. The default value is 100. //' @param N The number of points to sample. The default value is \eqn{5\cdot 10^5}. //' //' @references \cite{L. Cales, A. Chalkis, I.Z. Emiris, V. Fisikopoulos, @@ -36,19 +36,19 @@ //' h1 = h1 / 1000 //' h2=runif(n = 10, min = 1, max = 1000) //' h2 = h2 / 1000 -//' cop = copula(h1=h1, h2=h2, numSlices = 10, N = 100000) +//' cop = copula(R1 = h1, R2 = h2, M = 10, N = 100000) //' //' # compute a copula for a family of parallel hyperplanes and a family of conentric ellipsoids //' h = runif(n = 10, min = 1, max = 1000) //' h = h / 1000 //' E = replicate(10, rnorm(20)) //' E = cov(E) -//' cop = copula(h1=h, E=E, numSlices=10, N=100000) +//' cop = copula(R1 = h, Sigma = E, M = 10, N = 100000) //' //' @export // [[Rcpp::export]] -Rcpp::NumericMatrix copula (Rcpp::NumericVector h1 = R_NilValue, Rcpp::NumericVector h2 = R_NilValue, - Rcpp::NumericMatrix E = R_NilValue, Rcpp::Nullable numSlices = R_NilValue, +Rcpp::NumericMatrix copula (Rcpp::NumericVector R1 = R_NilValue, Rcpp::NumericVector R2 = R_NilValue, + Rcpp::NumericMatrix Sigma = R_NilValue, Rcpp::Nullable M = R_NilValue, Rcpp::Nullable N = R_NilValue){ typedef double NT; @@ -58,8 +58,8 @@ Rcpp::NumericMatrix copula (Rcpp::NumericVector h1 = R_NilValue, Rcpp::NumericVe typedef Eigen::Matrix MT; unsigned int num_slices = 100, numpoints = 500000; - if (numSlices.isNotNull()) { - num_slices = Rcpp::as(numSlices); + if (M.isNotNull()) { + num_slices = Rcpp::as(M); } if (N.isNotNull()) { @@ -70,22 +70,22 @@ Rcpp::NumericMatrix copula (Rcpp::NumericVector h1 = R_NilValue, Rcpp::NumericVe std::vector > StdCopula; unsigned int dim = Rcpp::as >(h1).size(), i, j; - if(!h1.isNotNull()) { + if(!R1.isNotNull()) { throw Rcpp::exception("You have to give at least one normal of a hyperplane!"); } - std::vector hyp1 = Rcpp::as >(h1); - if (h2.isNotNull()) { + std::vector hyp1 = Rcpp::as >(R1); + if (R2.isNotNull()) { - std::vector hyp2 = Rcpp::as < std::vector < NT > > (h2); + std::vector hyp2 = Rcpp::as < std::vector < NT > > (R2); StdCopula = twoParHypFam(dim, numpoints, num_slices, hyp1, hyp2); - } else if (E.isNotNull()) { + } else if (Sigma.isNotNull()) { std::vector > Gin(dim, std::vector(dim)); - MT EE = Rcpp::as(E); + MT EE = Rcpp::as(Sigma); for (i=0; i Date: Wed, 12 Feb 2020 21:52:34 +0200 Subject: [PATCH 05/13] fix bugs in copulas.cpp --- R-proj/R/RcppExports.R | 36 ++++--------- R-proj/R/compute_indicators.R | 6 +-- R-proj/src/copula.cpp | 21 +++++--- R-proj/src/copula2.cpp | 84 ------------------------------ include/convex_bodies/ellipsoids.h | 6 +-- 5 files changed, 29 insertions(+), 124 deletions(-) delete mode 100644 R-proj/src/copula2.cpp diff --git a/R-proj/R/RcppExports.R b/R-proj/R/RcppExports.R index 98c9cbf7f..b1001c82b 100644 --- a/R-proj/R/RcppExports.R +++ b/R-proj/R/RcppExports.R @@ -3,12 +3,13 @@ #' Construct a copula using uniform sampling from the unit simplex #' -#' Given two families of parallel hyperplanes intersecting the canonical simplex, this function uniformly samples from the canonical simplex and construct an approximation of the bivariate probability distribution, called copula. +#' Given two families of parallel hyperplanes or a family of parallel hyperplanes and a family of concentric ellispoids centered at the origin intersecting the canonical simplex, this function uniformly samples from the canonical simplex and construct an approximation of the bivariate probability distribution, called copula. #' #' @param h1 A \eqn{d}-dimensional vector that describes the direction of the first family of parallel hyperplanes. -#' @param h2 A \eqn{d}-dimensional vector that describes the direction of the second family of parallel hyperplanes. -#' @param numSlices The number of the slices for the copula. Default value is 100. -#' @param N The number of points to sample. Default value is \eqn{4\cdot 10^6}. +#' @param h2 Optional. A \eqn{d}-dimensional vector that describes the direction of the second family of parallel hyperplanes. +#' @param E Optional. The \eqn{d\times d} symmetric positive semidefine matrix that describes the family of concentric ellipsoids centered at the origin. +#' @param M The number of the slices for the copula. The default value is 100. +#' @param N The number of points to sample. The default value is \eqn{5\cdot 10^5}. #' #' @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.} @@ -20,35 +21,18 @@ #' h1 = h1 / 1000 #' h2=runif(n = 10, min = 1, max = 1000) #' h2 = h2 / 1000 -#' cop = copula1(h1=h1, h2=h2, numSlices = 10, N = 100000) -#' @export -copula1 <- function(h1, h2, numSlices, N) { - .Call(`_volesti_copula1`, h1, h2, numSlices, N) -} - -#' Construct a copula using uniform sampling from the unit simplex -#' -#' Given a family of parallel hyperplanes and a family of concentric ellispoids centered at the origin intersecting the canonical simplex, this function uniformly samples from the canonical simplex and construct an approximation of the bivariate probability distribution, called copula. -#' -#' @param h A \eqn{d}-dimensional vector that describes the direction of the first family of parallel hyperplanes. -#' @param E The \eqn{d\times d} symmetric positive semidefine matrix that describes the family of concentric ellipsoids centered at the origin. -#' @param numSlices The number of the slices for the copula. Default value is 100. -#' @param N The number of points to sample. Default value is \eqn{4\cdot 10^6}. +#' cop = copula(R1 = h1, R2 = h2, M = 10, N = 100000) #' -#' @references \cite{L. Cales, A. Chalkis, I.Z. Emiris, V. Fisikopoulos, -#' \dQuote{Practical volume computation of structured convex bodies, and an application to modeling portfolio dependencies and financial crises,} \emph{Proc. of Symposium on Computational Geometry, Budapest, Hungary,} 2018.} -#' -#' @return A \eqn{numSlices\times numSlices} numerical matrix that corresponds to a copula. -#' @examples #' # compute a copula for a family of parallel hyperplanes and a family of conentric ellipsoids #' h = runif(n = 10, min = 1, max = 1000) #' h = h / 1000 #' E = replicate(10, rnorm(20)) #' E = cov(E) -#' cop = copula2(h=h, E=E, numSlices=10, N=100000) +#' cop = copula(R1 = h, Sigma = E, M = 10, N = 100000) +#' #' @export -copula2 <- function(h, E, numSlices, N) { - .Call(`_volesti_copula2`, h, E, numSlices, N) +copula <- function(R1 = NULL, R2 = NULL, Sigma = NULL, M = NULL, N = NULL) { + .Call(`_volesti_copula`, R1, R2, Sigma, M, N) } #' Compute the exact volume of (a) a zonotope (b) an arbitrary simplex (c) a unit simplex (d) a cross polytope (e) a hypercube diff --git a/R-proj/R/compute_indicators.R b/R-proj/R/compute_indicators.R index eaac88c69..8fd7e8de5 100644 --- a/R-proj/R/compute_indicators.R +++ b/R-proj/R/compute_indicators.R @@ -26,7 +26,7 @@ compute_indicators <- function(MatReturns, W = NULL, M = NULL, N = NULL) { wl = W-1 indicators = c() - print(nrows-wl) + #print(nrows-wl) for (i in 1:(nrows-wl)) { Win=i:(i+wl) @@ -62,7 +62,7 @@ compute_indicators <- function(MatReturns, W = NULL, M = NULL, N = NULL) { } } indicators = c(indicators, blue_mass / red_mass) - print(length(indicators)) + #print(length(indicators)) } n = length(indicators) @@ -88,6 +88,6 @@ compute_indicators <- function(MatReturns, W = NULL, M = NULL, N = NULL) { } } - return(list("indicators" = indicatos, market_states = col) + return(list("indicators" = indicatos, market_states = col)) } diff --git a/R-proj/src/copula.cpp b/R-proj/src/copula.cpp index b79866d69..ee783435b 100644 --- a/R-proj/src/copula.cpp +++ b/R-proj/src/copula.cpp @@ -8,6 +8,7 @@ //Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 program. #include +#include #include #include "ellipsoids.h" #include "cartesian_geom/cartesian_kernel.h" @@ -47,15 +48,19 @@ //' //' @export // [[Rcpp::export]] -Rcpp::NumericMatrix copula (Rcpp::NumericVector R1 = R_NilValue, Rcpp::NumericVector R2 = R_NilValue, - Rcpp::NumericMatrix Sigma = R_NilValue, Rcpp::Nullable M = R_NilValue, - Rcpp::Nullable N = R_NilValue){ +Rcpp::NumericMatrix copula (Rcpp::Nullable R1 = R_NilValue, + Rcpp::Nullable R2 = R_NilValue, + Rcpp::Nullable Sigma = R_NilValue, + Rcpp::Nullable M = R_NilValue, + Rcpp::Nullable N = R_NilValue){ typedef double NT; typedef Cartesian Kernel; typedef typename Kernel::Point Point; typedef boost::mt19937 RNGType; typedef Eigen::Matrix MT; + typedef Eigen::Matrix VT; + typedef copula_ellipsoid CopEll; unsigned int num_slices = 100, numpoints = 500000; if (M.isNotNull()) { @@ -68,7 +73,7 @@ Rcpp::NumericMatrix copula (Rcpp::NumericVector R1 = R_NilValue, Rcpp::NumericVe Rcpp::NumericMatrix copula(num_slices, num_slices); std::vector > StdCopula; - unsigned int dim = Rcpp::as >(h1).size(), i, j; + unsigned int dim = Rcpp::as >(R1).size(), i, j; if(!R1.isNotNull()) { @@ -86,8 +91,8 @@ Rcpp::NumericMatrix copula (Rcpp::NumericVector R1 = R_NilValue, Rcpp::NumericVe std::vector > Gin(dim, std::vector(dim)); MT EE = Rcpp::as(Sigma); - for (i=0; i -#include -#include "ellipsoids.h" -#include -#include "cartesian_geom/cartesian_kernel.h" -#include -#include -#include "simplex_samplers.h" -#include "copulas.h" - -//' Construct a copula using uniform sampling from the unit simplex -//' -//' Given a family of parallel hyperplanes and a family of concentric ellispoids centered at the origin intersecting the canonical simplex, this function uniformly samples from the canonical simplex and construct an approximation of the bivariate probability distribution, called copula. -//' -//' @param h A \eqn{d}-dimensional vector that describes the direction of the first family of parallel hyperplanes. -//' @param E The \eqn{d\times d} symmetric positive semidefine matrix that describes the family of concentric ellipsoids centered at the origin. -//' @param numSlices The number of the slices for the copula. Default value is 100. -//' @param N The number of points to sample. Default value is \eqn{4\cdot 10^6}. -//' -//' @references \cite{L. Cales, A. Chalkis, I.Z. Emiris, V. Fisikopoulos, -//' \dQuote{Practical volume computation of structured convex bodies, and an application to modeling portfolio dependencies and financial crises,} \emph{Proc. of Symposium on Computational Geometry, Budapest, Hungary,} 2018.} -//' -//' @return A \eqn{numSlices\times numSlices} numerical matrix that corresponds to a copula. -//' @examples -//' # compute a copula for a family of parallel hyperplanes and a family of conentric ellipsoids -//' h = runif(n = 10, min = 1, max = 1000) -//' h = h / 1000 -//' E = replicate(10, rnorm(20)) -//' E = cov(E) -//' cop = copula2(h=h, E=E, numSlices=10, N=100000) -//' @export -// [[Rcpp::export]] -Rcpp::NumericMatrix copula2 (Rcpp::NumericVector h, Rcpp::NumericMatrix E, - Rcpp::Nullable numSlices, Rcpp::Nullable N){ - - typedef double NT; - typedef Cartesian Kernel; - typedef typename Kernel::Point Point; - typedef boost::mt19937 RNGType; - typedef copula_ellipsoid CopEll; - - unsigned int num_slices = 100, numpoints = 4000000; - - if (numSlices.isNotNull()) { - num_slices = Rcpp::as(numSlices); - } - - if (N.isNotNull()) { - numpoints = Rcpp::as(N); - } - - Rcpp::NumericMatrix copula(num_slices, num_slices); - std::vector > StdCopula; - unsigned int dim = h.size(), i, j; - - std::vector hyp1 = Rcpp::as >(h);; - std::vector > Gin(dim, std::vector(dim)); - - for (i=0; i(dim, numpoints, num_slices, hyp1, Ell); - - for(i=0; i -template +template class copula_ellipsoid{ private: typedef typename Point::FT NT; typedef typename std::vector::iterator viterator; - typedef Eigen::Matrix MT; - typedef Eigen::Matrix VT; + //typedef Eigen::Matrix MT; + //typedef Eigen::Matrix VT; MT G; unsigned int dim; public: From 8693492f3df30684533224de584b5a3247cecb00 Mon Sep 17 00:00:00 2001 From: tolischal Date: Thu, 13 Feb 2020 15:01:27 +0200 Subject: [PATCH 06/13] fix bugs in compute_indicators.R --- R-proj/NAMESPACE | 2 +- R-proj/R/compute_indicators.R | 16 ++++++------ R-proj/man/compute_indicators.Rd | 29 ++++++++++++++++++++ R-proj/man/copula.Rd | 45 ++++++++++++++++++++++++++++++++ R-proj/man/copula1.Rd | 35 ------------------------- R-proj/man/copula2.Rd | 35 ------------------------- 6 files changed, 83 insertions(+), 79 deletions(-) create mode 100644 R-proj/man/compute_indicators.Rd create mode 100644 R-proj/man/copula.Rd delete mode 100644 R-proj/man/copula1.Rd delete mode 100644 R-proj/man/copula2.Rd diff --git a/R-proj/NAMESPACE b/R-proj/NAMESPACE index 9710f5374..629a87d1d 100644 --- a/R-proj/NAMESPACE +++ b/R-proj/NAMESPACE @@ -1,7 +1,7 @@ # Generated by roxygen2: do not edit by hand -export(copula) export(compute_indicators) +export(copula) export(exact_vol) export(file_to_polytope) export(frustum_of_simplex) diff --git a/R-proj/R/compute_indicators.R b/R-proj/R/compute_indicators.R index 8fd7e8de5..3201bd534 100644 --- a/R-proj/R/compute_indicators.R +++ b/R-proj/R/compute_indicators.R @@ -30,11 +30,11 @@ compute_indicators <- function(MatReturns, W = NULL, M = NULL, N = NULL) { for (i in 1:(nrows-wl)) { Win=i:(i+wl) - if("tawny" %in% rownames(installed.packages()) == FALSE) { - E = cov(MatReturns[Win,]) - } else { - E = tawny::cov.shrink(MatReturns[Win,]) - } + #if("tawny" %in% rownames(installed.packages()) == FALSE && FALSE) { + E = cov(MatReturns[Win,]) + #} else { + #E = tawny::cov.shrink(MatReturns[Win,]) + #} compRet = rep(1,nassets) for (j in 1:nassets) { @@ -61,8 +61,8 @@ compute_indicators <- function(MatReturns, W = NULL, M = NULL, N = NULL) { } } } - indicators = c(indicators, blue_mass / red_mass) - #print(length(indicators)) + #indicators = c(indicators, blue_mass / red_mass) + print(length(indicators)) } n = length(indicators) @@ -88,6 +88,6 @@ compute_indicators <- function(MatReturns, W = NULL, M = NULL, N = NULL) { } } - return(list("indicators" = indicatos, market_states = col)) + return(list("indicators" = indicators, market_states = col)) } diff --git a/R-proj/man/compute_indicators.Rd b/R-proj/man/compute_indicators.Rd new file mode 100644 index 000000000..dbc267597 --- /dev/null +++ b/R-proj/man/compute_indicators.Rd @@ -0,0 +1,29 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/compute_indicators.R +\name{compute_indicators} +\alias{compute_indicators} +\title{Compute an indicator for each time period that describes the state of a market.} +\usage{ +compute_indicators(MatReturns, W = NULL, M = NULL, N = NULL) +} +\arguments{ +\item{MatReturns}{A \eqn{d}-dimensional vector that describes the direction of the first family of parallel hyperplanes.} + +\item{W}{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}.} +} +\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 W, this function computes an approximation of the joint distribution (copula) between portfolios' return and volatility in each time period implied by W. +For each copula it computes an indicator: large value corresponds to a crisis period and a small value to a normal period. +The periods over which the indicator is greater than 1 for more than 60 consecutives sliding windows are warnings and for more than 100 are crisis. The sliding window is shifted by one day. +} +\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.} +} diff --git a/R-proj/man/copula.Rd b/R-proj/man/copula.Rd new file mode 100644 index 000000000..f2bec2861 --- /dev/null +++ b/R-proj/man/copula.Rd @@ -0,0 +1,45 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{copula} +\alias{copula} +\title{Construct a copula using uniform sampling from the unit simplex} +\usage{ +copula(R1 = NULL, R2 = NULL, Sigma = NULL, M = NULL, N = NULL) +} +\arguments{ +\item{M}{The number of the 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{h1}{A \eqn{d}-dimensional vector that describes the direction of the first family of parallel hyperplanes.} + +\item{h2}{Optional. A \eqn{d}-dimensional vector that describes the direction of the second family of parallel hyperplanes.} + +\item{E}{Optional. The \eqn{d\times d} symmetric positive semidefine matrix that describes the family of concentric ellipsoids centered at the origin.} +} +\value{ +A \eqn{numSlices\times numSlices} numerical matrix that corresponds to a copula. +} +\description{ +Given two families of parallel hyperplanes or a family of parallel hyperplanes and a family of concentric ellispoids centered at the origin intersecting the canonical simplex, this function uniformly samples from the canonical simplex and construct an approximation of the bivariate probability distribution, called copula. +} +\examples{ +# compute a copula for two random families of parallel hyperplanes +h1 = runif(n = 10, min = 1, max = 1000) +h1 = h1 / 1000 +h2=runif(n = 10, min = 1, max = 1000) +h2 = h2 / 1000 +cop = copula(R1 = h1, R2 = h2, M = 10, N = 100000) + +# compute a copula for a family of parallel hyperplanes and a family of conentric ellipsoids +h = runif(n = 10, min = 1, max = 1000) +h = h / 1000 +E = replicate(10, rnorm(20)) +E = cov(E) +cop = copula(R1 = h, Sigma = E, M = 10, N = 100000) + +} +\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.} +} diff --git a/R-proj/man/copula1.Rd b/R-proj/man/copula1.Rd deleted file mode 100644 index a9413bc76..000000000 --- a/R-proj/man/copula1.Rd +++ /dev/null @@ -1,35 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/RcppExports.R -\name{copula1} -\alias{copula1} -\title{Construct a copula using uniform sampling from the unit simplex} -\usage{ -copula1(h1, h2, numSlices, N) -} -\arguments{ -\item{h1}{A \eqn{d}-dimensional vector that describes the direction of the first family of parallel hyperplanes.} - -\item{h2}{A \eqn{d}-dimensional vector that describes the direction of the second family of parallel hyperplanes.} - -\item{numSlices}{The number of the slices for the copula. Default value is 100.} - -\item{N}{The number of points to sample. Default value is \eqn{4\cdot 10^6}.} -} -\value{ -A \eqn{numSlices\times numSlices} numerical matrix that corresponds to a copula. -} -\description{ -Given two families of parallel hyperplanes intersecting the canonical simplex, this function uniformly samples from the canonical simplex and construct an approximation of the bivariate probability distribution, called copula. -} -\examples{ -# compute a copula for two random families of parallel hyperplanes -h1 = runif(n = 10, min = 1, max = 1000) -h1 = h1 / 1000 -h2=runif(n = 10, min = 1, max = 1000) -h2 = h2 / 1000 -cop = copula1(h1=h1, h2=h2, numSlices = 10, N = 100000) -} -\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.} -} diff --git a/R-proj/man/copula2.Rd b/R-proj/man/copula2.Rd deleted file mode 100644 index b693df821..000000000 --- a/R-proj/man/copula2.Rd +++ /dev/null @@ -1,35 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/RcppExports.R -\name{copula2} -\alias{copula2} -\title{Construct a copula using uniform sampling from the unit simplex} -\usage{ -copula2(h, E, numSlices, N) -} -\arguments{ -\item{h}{A \eqn{d}-dimensional vector that describes the direction of the first family of parallel hyperplanes.} - -\item{E}{The \eqn{d\times d} symmetric positive semidefine matrix that describes the family of concentric ellipsoids centered at the origin.} - -\item{numSlices}{The number of the slices for the copula. Default value is 100.} - -\item{N}{The number of points to sample. Default value is \eqn{4\cdot 10^6}.} -} -\value{ -A \eqn{numSlices\times numSlices} numerical matrix that corresponds to a copula. -} -\description{ -Given a family of parallel hyperplanes and a family of concentric ellispoids centered at the origin intersecting the canonical simplex, this function uniformly samples from the canonical simplex and construct an approximation of the bivariate probability distribution, called copula. -} -\examples{ -# compute a copula for a family of parallel hyperplanes and a family of conentric ellipsoids -h = runif(n = 10, min = 1, max = 1000) -h = h / 1000 -E = replicate(10, rnorm(20)) -E = cov(E) -cop = copula2(h=h, E=E, numSlices=10, N=100000) -} -\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.} -} From fc17489469b05cad9589b4b072f0830159506827 Mon Sep 17 00:00:00 2001 From: tolischal Date: Thu, 13 Feb 2020 15:41:04 +0200 Subject: [PATCH 07/13] fix R tests for copulas and set to lowercase the inputvariables in the new functions --- R-proj/DESCRIPTION | 4 +-- R-proj/NAMESPACE | 1 + R-proj/R/RcppExports.R | 18 +++++----- R-proj/R/compute_indicators.R | 50 ++++++++++++++-------------- R-proj/man/compute_indicators.Rd | 10 +++--- R-proj/man/copula.Rd | 16 ++++----- R-proj/src/copula.cpp | 46 ++++++++++++------------- R-proj/tests/testthat/test_copulas.R | 4 +-- 8 files changed, 75 insertions(+), 74 deletions(-) diff --git a/R-proj/DESCRIPTION b/R-proj/DESCRIPTION index d3188e83a..c0d95f7fd 100644 --- a/R-proj/DESCRIPTION +++ b/R-proj/DESCRIPTION @@ -15,9 +15,9 @@ Version: 1.1.0 Date: 2020-01-23 Maintainer: Vissarion Fisikopoulos Depends: Rcpp (>= 0.12.17) -Imports: methods +Imports: methods, tawny LinkingTo: Rcpp, RcppEigen, BH -Suggests: testthat, tawny +Suggests: testthat Encoding: UTF-8 RoxygenNote: 6.1.1 BugReports: https://github.com/GeomScale/volume_approximation/issues diff --git a/R-proj/NAMESPACE b/R-proj/NAMESPACE index 629a87d1d..f564cc008 100644 --- a/R-proj/NAMESPACE +++ b/R-proj/NAMESPACE @@ -21,6 +21,7 @@ export(volume) exportPattern("^[[:alpha:]]+") importFrom("methods","new") importFrom("utils","read.csv") +importFrom("tawny", "cov.shrink") importFrom(Rcpp,evalCpp) importFrom(Rcpp,loadModule) useDynLib(volesti, .registration=TRUE) diff --git a/R-proj/R/RcppExports.R b/R-proj/R/RcppExports.R index b1001c82b..d6704588e 100644 --- a/R-proj/R/RcppExports.R +++ b/R-proj/R/RcppExports.R @@ -5,11 +5,11 @@ #' #' Given two families of parallel hyperplanes or a family of parallel hyperplanes and a family of concentric ellispoids centered at the origin intersecting the canonical simplex, this function uniformly samples from the canonical simplex and construct an approximation of the bivariate probability distribution, called copula. #' -#' @param h1 A \eqn{d}-dimensional vector that describes the direction of the first family of parallel hyperplanes. -#' @param h2 Optional. A \eqn{d}-dimensional vector that describes the direction of the second family of parallel hyperplanes. -#' @param E Optional. The \eqn{d\times d} symmetric positive semidefine matrix that describes the family of concentric ellipsoids centered at the origin. -#' @param M The number of the slices for the copula. The default value is 100. -#' @param N The number of points to sample. The default value is \eqn{5\cdot 10^5}. +#' @param r1 A \eqn{d}-dimensional vector that describes the direction of the first family of parallel hyperplanes. +#' @param r2 Optional. A \eqn{d}-dimensional vector that describes the direction of the second family of parallel hyperplanes. +#' @param sigma Optional. The \eqn{d\times d} symmetric positive semidefine matrix that describes the family of concentric ellipsoids centered at the origin. +#' @param m The number of the slices for the copula. The default value is 100. +#' @param n The number of points to sample. The default value is \eqn{5\cdot 10^5}. #' #' @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.} @@ -21,18 +21,18 @@ #' h1 = h1 / 1000 #' h2=runif(n = 10, min = 1, max = 1000) #' h2 = h2 / 1000 -#' cop = copula(R1 = h1, R2 = h2, M = 10, N = 100000) +#' cop = copula(r1 = h1, r2 = h2, m = 10, n = 100000) #' #' # compute a copula for a family of parallel hyperplanes and a family of conentric ellipsoids #' h = runif(n = 10, min = 1, max = 1000) #' h = h / 1000 #' E = replicate(10, rnorm(20)) #' E = cov(E) -#' cop = copula(R1 = h, Sigma = E, M = 10, N = 100000) +#' cop = copula(r1 = h, sigma = E, m = 10, n = 100000) #' #' @export -copula <- function(R1 = NULL, R2 = NULL, Sigma = NULL, M = NULL, N = NULL) { - .Call(`_volesti_copula`, R1, R2, Sigma, M, N) +copula <- function(r1 = NULL, r2 = NULL, sigma = NULL, m = NULL, n = NULL) { + .Call(`_volesti_copula`, r1, r2, sigma, m, n) } #' Compute the exact volume of (a) a zonotope (b) an arbitrary simplex (c) a unit simplex (d) a cross polytope (e) a hypercube diff --git a/R-proj/R/compute_indicators.R b/R-proj/R/compute_indicators.R index 3201bd534..1203a6f7a 100644 --- a/R-proj/R/compute_indicators.R +++ b/R-proj/R/compute_indicators.R @@ -4,10 +4,10 @@ #' For each copula it computes an indicator: large value corresponds to a crisis period and a small value to a normal period. #' The periods over which the indicator is greater than 1 for more than 60 consecutives sliding windows are warnings and for more than 100 are crisis. The sliding window is shifted by one day. #' -#' @param MatReturns A \eqn{d}-dimensional vector that describes the direction of the first family of parallel hyperplanes. -#' @param W 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 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}. #' #' @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.} @@ -15,15 +15,15 @@ #' @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. #' #' @export -compute_indicators <- function(MatReturns, W = NULL, M = NULL, N = NULL) { +compute_indicators <- function(returns, win_length = NULL, m = NULL, n = NULL) { - if (is.null(W)) W = 60 - if (is.null(M)) M = 100 - if (is.null(N)) N = 500000 + if (is.null(win_length)) win_length = 60 + if (is.null(m)) m = 100 + if (is.null(n)) n = 500000 - nrows = dim(MatReturns)[1] - nassets = dim(MatReturns)[2] - wl = W-1 + nrows = dim(returns)[1] + nassets = dim(returns)[2] + wl = win_length-1 indicators = c() #print(nrows-wl) @@ -31,47 +31,47 @@ compute_indicators <- function(MatReturns, W = NULL, M = NULL, N = NULL) { Win=i:(i+wl) #if("tawny" %in% rownames(installed.packages()) == FALSE && FALSE) { - E = cov(MatReturns[Win,]) + #E = cov(returns[Win,]) #} else { - #E = tawny::cov.shrink(MatReturns[Win,]) + E = tawny::cov.shrink(returns[Win,]) #} compRet = rep(1,nassets) for (j in 1:nassets) { for (k in Win) { - compRet[j] = compRet[j] * (1 + MatReturns[k, j]) + compRet[j] = compRet[j] * (1 + returns[k, j]) } compRet[j] = compRet[j] - 1 } - cop = copula(R1 = compRet, Sigma = E, M = M, N = N) + cop = copula(r1 = compRet, sigma = E, m = m, n = n) 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) { - if (row+col<0.8*M || row+col>1.2*M) { + for (row in 1:m) { + for (col in 1:m) { + if (row-col<=0.2*m && row-col>=-0.2*m) { + if (row+col<0.8*m || row+col>1.2*m) { red_mass = red_mass + cop[row,col] } } else { - if (row+col>=0.8*M && row+col<=1.2*M) { + if (row+col>=0.8*m && row+col<=1.2*m) { blue_mass = blue_mass + cop[row,col] } } } } - #indicators = c(indicators, blue_mass / red_mass) - print(length(indicators)) + indicators = c(indicators, blue_mass / red_mass) + #print(length(indicators)) } - n = length(indicators) + N = length(indicators) index = 0 set_index = FALSE - col = rep("normal",n) + col = rep("normal", N) - for (i in 1:n) { + for (i in 1:N) { if(indicators[i]>1 && !set_index){ index = i diff --git a/R-proj/man/compute_indicators.Rd b/R-proj/man/compute_indicators.Rd index dbc267597..bf9d720b4 100644 --- a/R-proj/man/compute_indicators.Rd +++ b/R-proj/man/compute_indicators.Rd @@ -4,16 +4,16 @@ \alias{compute_indicators} \title{Compute an indicator for each time period that describes the state of a market.} \usage{ -compute_indicators(MatReturns, W = NULL, M = NULL, N = NULL) +compute_indicators(returns, win_length = NULL, m = NULL, n = NULL) } \arguments{ -\item{MatReturns}{A \eqn{d}-dimensional vector that describes the direction of the first family of parallel hyperplanes.} +\item{returns}{A \eqn{d}-dimensional vector that describes the direction of the first family of parallel hyperplanes.} -\item{W}{Optional. The length of the sliding window. The default value is 60.} +\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{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{n}{Optional. The number of points to sample. The default value is \eqn{5\cdot 10^5}.} } \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. diff --git a/R-proj/man/copula.Rd b/R-proj/man/copula.Rd index f2bec2861..5743fc616 100644 --- a/R-proj/man/copula.Rd +++ b/R-proj/man/copula.Rd @@ -4,18 +4,18 @@ \alias{copula} \title{Construct a copula using uniform sampling from the unit simplex} \usage{ -copula(R1 = NULL, R2 = NULL, Sigma = NULL, M = NULL, N = NULL) +copula(r1 = NULL, r2 = NULL, sigma = NULL, m = NULL, n = NULL) } \arguments{ -\item{M}{The number of the slices for the copula. The default value is 100.} +\item{r1}{A \eqn{d}-dimensional vector that describes the direction of the first family of parallel hyperplanes.} -\item{N}{The number of points to sample. The default value is \eqn{5\cdot 10^5}.} +\item{r2}{Optional. A \eqn{d}-dimensional vector that describes the direction of the second family of parallel hyperplanes.} -\item{h1}{A \eqn{d}-dimensional vector that describes the direction of the first family of parallel hyperplanes.} +\item{sigma}{Optional. The \eqn{d\times d} symmetric positive semidefine matrix that describes the family of concentric ellipsoids centered at the origin.} -\item{h2}{Optional. A \eqn{d}-dimensional vector that describes the direction of the second family of parallel hyperplanes.} +\item{m}{The number of the slices for the copula. The default value is 100.} -\item{E}{Optional. The \eqn{d\times d} symmetric positive semidefine matrix that describes the family of concentric ellipsoids centered at the origin.} +\item{n}{The number of points to sample. The default value is \eqn{5\cdot 10^5}.} } \value{ A \eqn{numSlices\times numSlices} numerical matrix that corresponds to a copula. @@ -29,14 +29,14 @@ h1 = runif(n = 10, min = 1, max = 1000) h1 = h1 / 1000 h2=runif(n = 10, min = 1, max = 1000) h2 = h2 / 1000 -cop = copula(R1 = h1, R2 = h2, M = 10, N = 100000) +cop = copula(r1 = h1, r2 = h2, m = 10, n = 100000) # compute a copula for a family of parallel hyperplanes and a family of conentric ellipsoids h = runif(n = 10, min = 1, max = 1000) h = h / 1000 E = replicate(10, rnorm(20)) E = cov(E) -cop = copula(R1 = h, Sigma = E, M = 10, N = 100000) +cop = copula(r1 = h, sigma = E, m = 10, n = 100000) } \references{ diff --git a/R-proj/src/copula.cpp b/R-proj/src/copula.cpp index ee783435b..8f3f7d805 100644 --- a/R-proj/src/copula.cpp +++ b/R-proj/src/copula.cpp @@ -21,11 +21,11 @@ //' //' Given two families of parallel hyperplanes or a family of parallel hyperplanes and a family of concentric ellispoids centered at the origin intersecting the canonical simplex, this function uniformly samples from the canonical simplex and construct an approximation of the bivariate probability distribution, called copula. //' -//' @param h1 A \eqn{d}-dimensional vector that describes the direction of the first family of parallel hyperplanes. -//' @param h2 Optional. A \eqn{d}-dimensional vector that describes the direction of the second family of parallel hyperplanes. -//' @param E Optional. The \eqn{d\times d} symmetric positive semidefine matrix that describes the family of concentric ellipsoids centered at the origin. -//' @param M The number of the slices for the copula. The default value is 100. -//' @param N The number of points to sample. The default value is \eqn{5\cdot 10^5}. +//' @param r1 A \eqn{d}-dimensional vector that describes the direction of the first family of parallel hyperplanes. +//' @param r2 Optional. A \eqn{d}-dimensional vector that describes the direction of the second family of parallel hyperplanes. +//' @param sigma Optional. The \eqn{d\times d} symmetric positive semidefine matrix that describes the family of concentric ellipsoids centered at the origin. +//' @param m The number of the slices for the copula. The default value is 100. +//' @param n The number of points to sample. The default value is \eqn{5\cdot 10^5}. //' //' @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.} @@ -37,22 +37,22 @@ //' h1 = h1 / 1000 //' h2=runif(n = 10, min = 1, max = 1000) //' h2 = h2 / 1000 -//' cop = copula(R1 = h1, R2 = h2, M = 10, N = 100000) +//' cop = copula(r1 = h1, r2 = h2, m = 10, n = 100000) //' //' # compute a copula for a family of parallel hyperplanes and a family of conentric ellipsoids //' h = runif(n = 10, min = 1, max = 1000) //' h = h / 1000 //' E = replicate(10, rnorm(20)) //' E = cov(E) -//' cop = copula(R1 = h, Sigma = E, M = 10, N = 100000) +//' cop = copula(r1 = h, sigma = E, m = 10, n = 100000) //' //' @export // [[Rcpp::export]] -Rcpp::NumericMatrix copula (Rcpp::Nullable R1 = R_NilValue, - Rcpp::Nullable R2 = R_NilValue, - Rcpp::Nullable Sigma = R_NilValue, - Rcpp::Nullable M = R_NilValue, - Rcpp::Nullable N = R_NilValue){ +Rcpp::NumericMatrix copula (Rcpp::Nullable r1 = R_NilValue, + Rcpp::Nullable r2 = R_NilValue, + Rcpp::Nullable sigma = R_NilValue, + Rcpp::Nullable m = R_NilValue, + Rcpp::Nullable n = R_NilValue){ typedef double NT; typedef Cartesian Kernel; @@ -63,34 +63,34 @@ Rcpp::NumericMatrix copula (Rcpp::Nullable R1 = R_NilValue, typedef copula_ellipsoid CopEll; unsigned int num_slices = 100, numpoints = 500000; - if (M.isNotNull()) { - num_slices = Rcpp::as(M); + if (m.isNotNull()) { + num_slices = Rcpp::as(m); } - if (N.isNotNull()) { - numpoints = Rcpp::as(N); + if (n.isNotNull()) { + numpoints = Rcpp::as(n); } Rcpp::NumericMatrix copula(num_slices, num_slices); std::vector > StdCopula; - unsigned int dim = Rcpp::as >(R1).size(), i, j; + unsigned int dim = Rcpp::as >(r1).size(), i, j; - if(!R1.isNotNull()) { + if(!r1.isNotNull()) { throw Rcpp::exception("You have to give at least one normal of a hyperplane!"); } - std::vector hyp1 = Rcpp::as >(R1); - if (R2.isNotNull()) { + std::vector hyp1 = Rcpp::as >(r1); + if (r2.isNotNull()) { - std::vector hyp2 = Rcpp::as < std::vector < NT > > (R2); + std::vector hyp2 = Rcpp::as < std::vector < NT > > (r2); StdCopula = twoParHypFam(dim, numpoints, num_slices, hyp1, hyp2); - } else if (Sigma.isNotNull()) { + } else if (sigma.isNotNull()) { std::vector > Gin(dim, std::vector(dim)); - MT EE = Rcpp::as(Sigma); + MT EE = Rcpp::as(sigma); for (int i=0; i Date: Fri, 14 Feb 2020 13:15:31 +0200 Subject: [PATCH 08/13] improve compute_indicators.R --- R-proj/DESCRIPTION | 2 +- R-proj/NAMESPACE | 2 +- R-proj/R/compute_indicators.R | 10 ++-------- 3 files changed, 4 insertions(+), 10 deletions(-) diff --git a/R-proj/DESCRIPTION b/R-proj/DESCRIPTION index c0d95f7fd..ab21adf38 100644 --- a/R-proj/DESCRIPTION +++ b/R-proj/DESCRIPTION @@ -15,7 +15,7 @@ Version: 1.1.0 Date: 2020-01-23 Maintainer: Vissarion Fisikopoulos Depends: Rcpp (>= 0.12.17) -Imports: methods, tawny +Imports: methods LinkingTo: Rcpp, RcppEigen, BH Suggests: testthat Encoding: UTF-8 diff --git a/R-proj/NAMESPACE b/R-proj/NAMESPACE index f564cc008..181f091f3 100644 --- a/R-proj/NAMESPACE +++ b/R-proj/NAMESPACE @@ -21,7 +21,7 @@ export(volume) exportPattern("^[[:alpha:]]+") importFrom("methods","new") importFrom("utils","read.csv") -importFrom("tawny", "cov.shrink") +importFrom("stats", "cov") importFrom(Rcpp,evalCpp) importFrom(Rcpp,loadModule) useDynLib(volesti, .registration=TRUE) diff --git a/R-proj/R/compute_indicators.R b/R-proj/R/compute_indicators.R index 1203a6f7a..fd859ec18 100644 --- a/R-proj/R/compute_indicators.R +++ b/R-proj/R/compute_indicators.R @@ -26,15 +26,10 @@ compute_indicators <- function(returns, win_length = NULL, m = NULL, n = NULL) { wl = win_length-1 indicators = c() - #print(nrows-wl) for (i in 1:(nrows-wl)) { Win=i:(i+wl) - #if("tawny" %in% rownames(installed.packages()) == FALSE && FALSE) { - #E = cov(returns[Win,]) - #} else { - E = tawny::cov.shrink(returns[Win,]) - #} + E = cov(returns[Win,]) compRet = rep(1,nassets) for (j in 1:nassets) { @@ -55,14 +50,13 @@ compute_indicators <- function(returns, win_length = NULL, m = NULL, n = NULL) { red_mass = red_mass + cop[row,col] } } else { - if (row+col>=0.8*m && row+col<=1.2*m) { + if (row+col>=0.8*m+1 && row+col<=1.2*m+1) { blue_mass = blue_mass + cop[row,col] } } } } indicators = c(indicators, blue_mass / red_mass) - #print(length(indicators)) } N = length(indicators) From 410f2ef1ce25a67ca6ec5d6e2dfc2dbdb744e34d Mon Sep 17 00:00:00 2001 From: Tolis Date: Fri, 14 Feb 2020 13:53:52 +0200 Subject: [PATCH 09/13] add volume field in R module polytope classes --- R-proj/src/polytopes_modules.cpp | 36 ++++++++++++++++++++++++++++---- 1 file changed, 32 insertions(+), 4 deletions(-) diff --git a/R-proj/src/polytopes_modules.cpp b/R-proj/src/polytopes_modules.cpp index 892138e98..753a472f0 100644 --- a/R-proj/src/polytopes_modules.cpp +++ b/R-proj/src/polytopes_modules.cpp @@ -12,8 +12,12 @@ class Hpolytope { Hpolytope(Rcpp::NumericMatrix _A, Rcpp::NumericVector _b) : A(_A), b(_b) { dimension = _A.ncol(); } + Hpolytope(Rcpp::NumericMatrix _A, Rcpp::NumericVector _b, double volume) : A(_A), b(_b), vol(volume) { + dimension = _A.ncol(); + } int type = 1; unsigned int dimension; + double vol = -1.0; Rcpp::NumericMatrix A; Rcpp::NumericVector b; @@ -25,8 +29,12 @@ class Vpolytope { Vpolytope(Rcpp::NumericMatrix _V) : V(_V) { dimension = _V.ncol(); } + Vpolytope(Rcpp::NumericMatrix _V, double volume) : V(_V), vol(volume) { + dimension = _V.ncol(); + } int type = 2; unsigned int dimension; + double vol = -1.0; Rcpp::NumericMatrix V; }; @@ -37,8 +45,12 @@ class Zonotope { Zonotope(Rcpp::NumericMatrix _G) : G(_G) { dimension = _G.ncol(); } + Zonotope(Rcpp::NumericMatrix _G, double volume) : G(_G), vol(volume) { + dimension = _G.ncol(); + } int type = 3; unsigned int dimension; + double vol = -1.0; Rcpp::NumericMatrix G; }; @@ -49,8 +61,12 @@ class VPinterVP { VPinterVP(Rcpp::NumericMatrix _V1, Rcpp::NumericMatrix _V2) : V1(_V1), V2(_V2) { dimension = _V1.ncol(); } + VPinterVP(Rcpp::NumericMatrix _V1, Rcpp::NumericMatrix _V2, double volume) : V1(_V1), V2(_V2),vol(volume) { + dimension = _V1.ncol(); + } int type = 4; unsigned int dimension; + double vol = -1.0; Rcpp::NumericMatrix V1; Rcpp::NumericMatrix V2; @@ -68,7 +84,8 @@ RCPP_MODULE(yada){ //' @field A \eqn{m\times d} numerical matrix A //' @field b \eqn{m}-dimensional vector b //' @field type An integer that declares the representation of the polytope. For H-representation the default value is 1. It has not be given to the constructor. - //' @field An integer that declares the dimension of the polytope. It has not be given to the constructor. + //' @field dimension An integer that declares the dimension of the polytope. It has not be given to the constructor. + //' @field volume The volume of the polytope. //' //' @example //' # create a 2-d unit simplex @@ -80,9 +97,11 @@ RCPP_MODULE(yada){ // expose the default constructor .constructor() .constructor() + .constructor() .field( "type", &Hpolytope::type ) .field( "dimension", &Hpolytope::dimension ) + .field( "volume", &Hpolytope::vol ) .field( "b", &Hpolytope::b ) .field( "A", &Hpolytope::A ); @@ -92,7 +111,8 @@ RCPP_MODULE(yada){ //' //' @field V \eqn{m\times d} numerical matrix that contains the vertices row-wise //' @field type An integer that declares the representation of the polytope. For H-representation the default value is 2. It has not be given to the constructor. - //' @field An integer that declares the dimension of the polytope. It has not be given to the constructor. + //' @field dimension An integer that declares the dimension of the polytope. It has not be given to the constructor. + //' @field volume The volume of the polytope. //' //' @example //' # Create a 2-d cube in V-representation @@ -103,9 +123,11 @@ RCPP_MODULE(yada){ // expose the default constructor .constructor() .constructor() + .constructor() .field( "type", &Vpolytope::type ) .field( "dimension", &Vpolytope::dimension ) + .field( "volume", &Vpolytope::vol ) .field( "V", &Vpolytope::V ); //' An exposed C++ class to represent a zonotope @@ -114,7 +136,8 @@ RCPP_MODULE(yada){ //' //' @field G \eqn{m\times d} numerical matrix that contains the segments (or generators) row-wise //' @field type An integer that declares the representation of the polytope. For H-representation the default value is 3. It has not be given to the constructor. - //' @field An integer that declares the dimension of the polytope. It has not be given to the constructor. + //' @field dimension An integer that declares the dimension of the polytope. It has not be given to the constructor. + //' @field volume The volume of the polytope. //' //' @example //' # Create a 2-d zonotope with 4 generators @@ -125,9 +148,11 @@ RCPP_MODULE(yada){ // expose the default constructor .constructor() .constructor() + .constructor() .field( "type", &Zonotope::type ) .field( "dimension", &Zonotope::dimension ) + .field( "volume", &Zonotope::vol ) .field( "G", &Zonotope::G ); //' An exposed C++ class to represent an intersection of two V-polytopes @@ -136,7 +161,8 @@ RCPP_MODULE(yada){ //' //' @field G \eqn{m\times d} numerical matrix that contains the segments (or generators) row-wise //' @field type An integer that declares the representation of the polytope. For H-representation the default value is 4. It has not be given to the constructor. - //' @field An integer that declares the dimension of the polytope. It has not be given to the constructor. + //' @field dimension An integer that declares the dimension of the polytope. It has not be given to the constructor. + //' @field volume The volume of the polytope. //' //' @example //' # define the intwrsection of a 2-d simplex with a 2-d cross polytope @@ -148,9 +174,11 @@ RCPP_MODULE(yada){ // expose the default constructor .constructor() .constructor() + .constructor() .field( "type", &VPinterVP::type ) .field( "dimension", &VPinterVP::dimension ) + .field( "volume", &VPinterVP::vol ) .field( "V1", &VPinterVP::V1 ) .field( "V2", &VPinterVP::V2 ); } From 774a0db612f8c5bdd859a62fab4feb2eeb8ad63c Mon Sep 17 00:00:00 2001 From: Tolis Date: Fri, 14 Feb 2020 14:10:56 +0200 Subject: [PATCH 10/13] add volume declaration for known polytopes in R interface --- R-proj/R/gen_cross.R | 4 ++-- R-proj/R/gen_cube.R | 4 ++-- R-proj/R/gen_prod_simplex.R | 2 +- R-proj/R/gen_simplex.R | 4 ++-- R-proj/R/gen_skinny_cube.R | 2 +- 5 files changed, 8 insertions(+), 8 deletions(-) diff --git a/R-proj/R/gen_cross.R b/R-proj/R/gen_cross.R index cec9dde9a..026735ae0 100644 --- a/R-proj/R/gen_cross.R +++ b/R-proj/R/gen_cross.R @@ -32,9 +32,9 @@ gen_cross <- function(dimension, repr) { Mat = Mat[,-c(1)] if (Vpoly_gen) { - P = Vpolytope$new(Mat) + P = Vpolytope$new(Mat, 2^dimension / prod(1:dimension)) } else { - P = Hpolytope$new(-Mat, b) + P = Hpolytope$new(-Mat, b, 2^dimension / prod(1:dimension)) } return(P) diff --git a/R-proj/R/gen_cube.R b/R-proj/R/gen_cube.R index 69c17bd8f..cce62489c 100644 --- a/R-proj/R/gen_cube.R +++ b/R-proj/R/gen_cube.R @@ -31,9 +31,9 @@ gen_cube <- function(dimension, repr) { b = Mat[,1] Mat = Mat[,-c(1)] if (Vpoly_gen) { - P = Vpolytope$new(Mat) + P = Vpolytope$new(Mat, 2^dimension) } else { - P = Hpolytope$new(-Mat, b) + P = Hpolytope$new(-Mat, b, 2^dimension) } return(P) diff --git a/R-proj/R/gen_prod_simplex.R b/R-proj/R/gen_prod_simplex.R index 934c3f01c..95e506948 100644 --- a/R-proj/R/gen_prod_simplex.R +++ b/R-proj/R/gen_prod_simplex.R @@ -22,7 +22,7 @@ gen_prod_simplex <- function(dimension) { b = Mat[,1] Mat = Mat[,-c(1)] - P = Hpolytope$new(-Mat, b) + P = Hpolytope$new(-Mat, b, (1/prod(1:dimension))^2) return(P) diff --git a/R-proj/R/gen_simplex.R b/R-proj/R/gen_simplex.R index b0434df2d..b1f04d917 100644 --- a/R-proj/R/gen_simplex.R +++ b/R-proj/R/gen_simplex.R @@ -32,9 +32,9 @@ gen_simplex <- function(dimension, repr) { Mat = Mat[,-c(1)] if (Vpoly_gen) { - P = Vpolytope$new(Mat) + P = Vpolytope$new(Mat, 1/prod(1:dimension)) } else { - P = Hpolytope$new(-Mat, b) + P = Hpolytope$new(-Mat, b, 1/prod(1:dimension)) } return(P) diff --git a/R-proj/R/gen_skinny_cube.R b/R-proj/R/gen_skinny_cube.R index 60b88df05..2c2dbc67f 100644 --- a/R-proj/R/gen_skinny_cube.R +++ b/R-proj/R/gen_skinny_cube.R @@ -22,7 +22,7 @@ gen_skinny_cube <- function(dimension) { b = Mat[,1] Mat = Mat[,-c(1)] - P = Hpolytope$new(-Mat, b) + P = Hpolytope$new(-Mat, b, 2^(dimension -1)*200) return(P) } From d363040059b25c4f5955a699e00c072a5fef7fa5 Mon Sep 17 00:00:00 2001 From: Tolis Date: Fri, 14 Feb 2020 15:16:07 +0200 Subject: [PATCH 11/13] modify exact_vol.cpp in R interface --- R-proj/DESCRIPTION | 2 +- R-proj/R/RcppExports.R | 20 +++--- R-proj/man/exact_vol.Rd | 20 ++---- R-proj/src/exact_vol.cpp | 128 +++++++++++++++------------------------ 4 files changed, 62 insertions(+), 108 deletions(-) diff --git a/R-proj/DESCRIPTION b/R-proj/DESCRIPTION index ab21adf38..ffda1a9fa 100644 --- a/R-proj/DESCRIPTION +++ b/R-proj/DESCRIPTION @@ -15,7 +15,7 @@ Version: 1.1.0 Date: 2020-01-23 Maintainer: Vissarion Fisikopoulos Depends: Rcpp (>= 0.12.17) -Imports: methods +Imports: methods, stats LinkingTo: Rcpp, RcppEigen, BH Suggests: testthat Encoding: UTF-8 diff --git a/R-proj/R/RcppExports.R b/R-proj/R/RcppExports.R index d6704588e..58db215af 100644 --- a/R-proj/R/RcppExports.R +++ b/R-proj/R/RcppExports.R @@ -35,21 +35,14 @@ copula <- function(r1 = NULL, r2 = NULL, sigma = NULL, m = NULL, n = NULL) { .Call(`_volesti_copula`, r1, r2, sigma, m, n) } -#' Compute the exact volume of (a) a zonotope (b) an arbitrary simplex (c) a unit simplex (d) a cross polytope (e) a hypercube +#' Compute the exact volume of (a) a zonotope (b) an arbitrary simplex in V-representation or (c) if the volume is known and declared by the input object. #' #' Given a zonotope (as an object of class Zonotope), this function computes the sum of the absolute values of the determinants of all the \eqn{d \times d} submatrices of the \eqn{m\times d} matrix \eqn{G} that contains row-wise the \eqn{m} \eqn{d}-dimensional segments that define the zonotope. #' For an arbitrary simplex that is given in V-representation this function computes the absolute value of the determinant formed by the simplex's points assuming it is shifted to the origin. -#' For a \eqn{d}-dimensional unit simplex, hypercube or cross polytope this function computes the exact well known formulas. #' -#' @param P A zonotope or a simplex in V-representation. -#' @param body A string that declares the type of the body for the exact sampling: a) \code{'simplex'} for the unit simplex, b) \code{'cross'} for the cross polytope, c) \code{'hypersphere'} for the hypersphere, d) \code{'cube'} for the unit cube. -#' @param Parameters A list for the parameters of the methods: -#' \itemize{ -#' \item{\code{dimension} }{ An integer that declares the dimension when exact sampling is enabled for a simplex or a hypersphere.} -#' \item{\code{radius} }{ The radius of the \eqn{d}-dimensional hypersphere. Default value is \eqn{1}.} -#' } +#' @param P A polytope #' -#' @return The exact volume of the zonotope +#' @return The exact volume of the input polytope, for zonotopes, simplices in V-representation and polytopes with known exact volume #' @examples #' #' # compute the exact volume of a 5-dimensional zonotope defined by the Minkowski sum of 10 segments @@ -63,10 +56,11 @@ copula <- function(r1 = NULL, r2 = NULL, sigma = NULL, m = NULL, n = NULL) { #' } #' #' # compute the exact volume the 10-dimensional cross polytope -#' vol = exact_vol(body = "cross", Parameters = list("dimension" = 10)) +#' P = gen_cross(10,'V') +#' vol = exact_vol(P) #' @export -exact_vol <- function(P = NULL, body = NULL, Parameters = NULL) { - .Call(`_volesti_exact_vol`, P, body, Parameters) +exact_vol <- function(P) { + .Call(`_volesti_exact_vol`, P) } #' Compute the percentage of the volume of the unit simplex that is contained in the intersection of a half-space and the unit simplex. diff --git a/R-proj/man/exact_vol.Rd b/R-proj/man/exact_vol.Rd index cc6a430b8..b6a6a4899 100644 --- a/R-proj/man/exact_vol.Rd +++ b/R-proj/man/exact_vol.Rd @@ -2,28 +2,19 @@ % Please edit documentation in R/RcppExports.R \name{exact_vol} \alias{exact_vol} -\title{Compute the exact volume of (a) a zonotope (b) an arbitrary simplex (c) a unit simplex (d) a cross polytope (e) a hypercube} +\title{Compute the exact volume of (a) a zonotope (b) an arbitrary simplex in V-representation or (c) if the volume is known and declared by the input object.} \usage{ -exact_vol(P = NULL, body = NULL, Parameters = NULL) +exact_vol(P) } \arguments{ -\item{P}{A zonotope or a simplex in V-representation.} - -\item{body}{A string that declares the type of the body for the exact sampling: a) \code{'simplex'} for the unit simplex, b) \code{'cross'} for the cross polytope, c) \code{'hypersphere'} for the hypersphere, d) \code{'cube'} for the unit cube.} - -\item{Parameters}{A list for the parameters of the methods: -\itemize{ -\item{\code{dimension} }{ An integer that declares the dimension when exact sampling is enabled for a simplex or a hypersphere.} -\item{\code{radius} }{ The radius of the \eqn{d}-dimensional hypersphere. Default value is \eqn{1}.} -}} +\item{P}{A polytope} } \value{ -The exact volume of the zonotope +The exact volume of the input polytope, for zonotopes, simplices in V-representation and polytopes with known exact volume } \description{ Given a zonotope (as an object of class Zonotope), this function computes the sum of the absolute values of the determinants of all the \eqn{d \times d} submatrices of the \eqn{m\times d} matrix \eqn{G} that contains row-wise the \eqn{m} \eqn{d}-dimensional segments that define the zonotope. For an arbitrary simplex that is given in V-representation this function computes the absolute value of the determinant formed by the simplex's points assuming it is shifted to the origin. -For a \eqn{d}-dimensional unit simplex, hypercube or cross polytope this function computes the exact well known formulas. } \examples{ @@ -38,5 +29,6 @@ vol = exact_vol(P) } # compute the exact volume the 10-dimensional cross polytope -vol = exact_vol(body = "cross", Parameters = list("dimension" = 10)) +P = gen_cross(10,'V') +vol = exact_vol(P) } diff --git a/R-proj/src/exact_vol.cpp b/R-proj/src/exact_vol.cpp index 23c0a8f13..8a6121173 100644 --- a/R-proj/src/exact_vol.cpp +++ b/R-proj/src/exact_vol.cpp @@ -24,21 +24,14 @@ FT factorial(FT n) return (n == 1 || n == 0) ? 1 : factorial(n - 1) * n; } -//' Compute the exact volume of (a) a zonotope (b) an arbitrary simplex (c) a unit simplex (d) a cross polytope (e) a hypercube +//' Compute the exact volume of (a) a zonotope (b) an arbitrary simplex in V-representation or (c) if the volume is known and declared by the input object. //' //' Given a zonotope (as an object of class Zonotope), this function computes the sum of the absolute values of the determinants of all the \eqn{d \times d} submatrices of the \eqn{m\times d} matrix \eqn{G} that contains row-wise the \eqn{m} \eqn{d}-dimensional segments that define the zonotope. //' For an arbitrary simplex that is given in V-representation this function computes the absolute value of the determinant formed by the simplex's points assuming it is shifted to the origin. -//' For a \eqn{d}-dimensional unit simplex, hypercube or cross polytope this function computes the exact well known formulas. //' -//' @param P A zonotope or a simplex in V-representation. -//' @param body A string that declares the type of the body for the exact sampling: a) \code{'simplex'} for the unit simplex, b) \code{'cross'} for the cross polytope, c) \code{'hypersphere'} for the hypersphere, d) \code{'cube'} for the unit cube. -//' @param Parameters A list for the parameters of the methods: -//' \itemize{ -//' \item{\code{dimension} }{ An integer that declares the dimension when exact sampling is enabled for a simplex or a hypersphere.} -//' \item{\code{radius} }{ The radius of the \eqn{d}-dimensional hypersphere. Default value is \eqn{1}.} -//' } +//' @param P A polytope //' -//' @return The exact volume of the zonotope +//' @return The exact volume of the input polytope, for zonotopes, simplices in V-representation and polytopes with known exact volume //' @examples //' //' # compute the exact volume of a 5-dimensional zonotope defined by the Minkowski sum of 10 segments @@ -52,84 +45,59 @@ FT factorial(FT n) //' } //' //' # compute the exact volume the 10-dimensional cross polytope -//' vol = exact_vol(body = "cross", Parameters = list("dimension" = 10)) +//' P = gen_cross(10,'V') +//' vol = exact_vol(P) //' @export // [[Rcpp::export]] -double exact_vol(Rcpp::Nullable P = R_NilValue, Rcpp::Nullable body = R_NilValue, - Rcpp::Nullable Parameters = R_NilValue){ +double exact_vol(Rcpp::Nullable P) { typedef double NT; - typedef Cartesian Kernel; - typedef typename Kernel::Point Point; - typedef boost::mt19937 RNGType; - typedef Eigen::Matrix VT; - typedef Eigen::Matrix MT; - - typedef Zonotope zonotope; - typedef VPolytope Vpolytope; - zonotope ZP; - Vpolytope VP; - int type, dim; - NT vol, rad; - - if (P.isNotNull()) { - type = Rcpp::as(P).field("type"); - if (!body.isNotNull()) { - dim = Rcpp::as(P).field("dimension"); - if (type == 3) { - ZP.init(dim, Rcpp::as(Rcpp::as(P).field("G")), - VT::Ones(Rcpp::as(Rcpp::as(P).field("G")).rows())); - vol = exact_zonotope_vol(ZP); - } else if (type == 2) { - if (Rcpp::as(Rcpp::as(P).field("V")).rows() == Rcpp::as(Rcpp::as(P).field("V")).cols()+1) { - MT V = Rcpp::as(Rcpp::as(P).field("V")).transpose(); - VT v0 = V.col(dim); - VT V2 = V.block(0, 0, dim, dim); - V2 = V2.colwise() - v0; - - vol = std::abs(V2.determinant()); - vol = vol / factorial(NT(dim)); - } else { - throw Rcpp::exception("Not a simplex!"); - } - } else { - throw Rcpp::exception("Wrong input!"); - } - } else { - throw Rcpp::exception("When a polytope is given, input for body has to be null."); - } - } else { - if (body.isNotNull()) { - if (Parameters.isNotNull()) { - if (Rcpp::as(Parameters).containsElementNamed("dimension")) { - dim = Rcpp::as(Rcpp::as(Parameters)["dimension"]); - } else { - throw Rcpp::exception("You have to declare the dimension!"); - } - } else { - throw Rcpp::exception("You have to declare the dimension!"); - } - if (Rcpp::as(body).compare(std::string("simplex"))==0) { - vol = 1.0 / factorial(NT(dim)); - } else if (Rcpp::as(body).compare(std::string("cross"))==0) { - vol = std::pow(2.0, NT(dim)); - vol = vol / factorial(NT(dim)); - } else if (Rcpp::as(body).compare(std::string("cube"))==0) { - vol = std::pow(2.0, NT(dim)); - } else if (Rcpp::as(body).compare(std::string("hypersphere"))==0) { - if (Rcpp::as(Parameters).containsElementNamed("radius")) { - rad = Rcpp::as(Rcpp::as(Parameters)["radius"]); - } else { - throw Rcpp::exception("You have to declare the dimension!"); - } - vol = (std::pow(M_PI,dim/2.0)*(std::pow(rad, dim) ) ) / (tgamma(dim/2.0+1)); - } else { - throw Rcpp::exception("Unknown body!"); + typedef Cartesian Kernel; + typedef typename Kernel::Point Point; + typedef boost::mt19937 RNGType; + typedef Eigen::Matrix VT; + typedef Eigen::Matrix MT; + + if (NT(Rcpp::as(P).field("volume")) > 0.0) { + return NT(Rcpp::as(P).field("volume")); + } + + int type = Rcpp::as(P).field("type"), dim; + NT vol; + + if (type == 2) { + + typedef VPolytope Vpolytope; + Vpolytope VP; + dim = Rcpp::as(P).field("dimension"); + + if (Rcpp::as(Rcpp::as(P).field("V")).rows() == + Rcpp::as(Rcpp::as(P).field("V")).cols() + 1) { + + MT V = Rcpp::as(Rcpp::as(P).field("V")).transpose(), V2(dim,dim); + VT v0 = V.col(dim); + + for (int i = 0; i < dim; ++i) { + V2.col(i) = V.col(i) - v0; } + vol = std::abs(V2.determinant()) / factorial(NT(dim)); } else { - throw Rcpp::exception("When a polytope is null, input for specific body required."); + throw Rcpp::exception("Volume unknown!"); } + + } else if (type == 3) { + + typedef Zonotope zonotope; + zonotope ZP; + dim = Rcpp::as(P).field("dimension"); + + ZP.init(dim, Rcpp::as(Rcpp::as(P).field("G")), + VT::Ones(Rcpp::as(Rcpp::as(P).field("G")).rows())); + vol = exact_zonotope_vol(ZP); + + } else { + throw Rcpp::exception("Volume unknown!"); } return vol; From 77689f2108a60da8753e2a2c4938fa20cc2c8da5 Mon Sep 17 00:00:00 2001 From: Tolis Date: Fri, 14 Feb 2020 15:41:37 +0200 Subject: [PATCH 12/13] set NaN the volume initialization in Rcpp modules of polytopes, update NEWS.md file --- R-proj/NAMESPACE | 2 +- R-proj/src/exact_vol.cpp | 2 +- R-proj/src/polytopes_modules.cpp | 8 ++++---- cran_gen/NEWS.md | 2 ++ 4 files changed, 8 insertions(+), 6 deletions(-) diff --git a/R-proj/NAMESPACE b/R-proj/NAMESPACE index 181f091f3..317a83a2c 100644 --- a/R-proj/NAMESPACE +++ b/R-proj/NAMESPACE @@ -20,8 +20,8 @@ export(sample_points) export(volume) exportPattern("^[[:alpha:]]+") importFrom("methods","new") -importFrom("utils","read.csv") importFrom("stats", "cov") +importFrom("utils","read.csv") importFrom(Rcpp,evalCpp) importFrom(Rcpp,loadModule) useDynLib(volesti, .registration=TRUE) diff --git a/R-proj/src/exact_vol.cpp b/R-proj/src/exact_vol.cpp index 8a6121173..eb58a9840 100644 --- a/R-proj/src/exact_vol.cpp +++ b/R-proj/src/exact_vol.cpp @@ -76,7 +76,7 @@ double exact_vol(Rcpp::Nullable P) { MT V = Rcpp::as(Rcpp::as(P).field("V")).transpose(), V2(dim,dim); VT v0 = V.col(dim); - + for (int i = 0; i < dim; ++i) { V2.col(i) = V.col(i) - v0; } diff --git a/R-proj/src/polytopes_modules.cpp b/R-proj/src/polytopes_modules.cpp index 753a472f0..464d476e4 100644 --- a/R-proj/src/polytopes_modules.cpp +++ b/R-proj/src/polytopes_modules.cpp @@ -17,7 +17,7 @@ class Hpolytope { } int type = 1; unsigned int dimension; - double vol = -1.0; + double vol = std::numeric_limits::signaling_NaN(); Rcpp::NumericMatrix A; Rcpp::NumericVector b; @@ -34,7 +34,7 @@ class Vpolytope { } int type = 2; unsigned int dimension; - double vol = -1.0; + double vol = std::numeric_limits::signaling_NaN(); Rcpp::NumericMatrix V; }; @@ -50,7 +50,7 @@ class Zonotope { } int type = 3; unsigned int dimension; - double vol = -1.0; + double vol = std::numeric_limits::signaling_NaN(); Rcpp::NumericMatrix G; }; @@ -66,7 +66,7 @@ class VPinterVP { } int type = 4; unsigned int dimension; - double vol = -1.0; + double vol = std::numeric_limits::signaling_NaN(); Rcpp::NumericMatrix V1; Rcpp::NumericMatrix V2; diff --git a/cran_gen/NEWS.md b/cran_gen/NEWS.md index 583961222..2a5eaafc1 100644 --- a/cran_gen/NEWS.md +++ b/cran_gen/NEWS.md @@ -18,4 +18,6 @@ * New volume computation algorithm. * Billiard walk for uniform sampling. +* Modified exact volume computation function. +* Improved functionality for finance applications. * Improved names for functions and input variables. From e44ebb3f6a5d10fbad5e7e1e849a05185ee1fe48 Mon Sep 17 00:00:00 2001 From: Tolis Date: Fri, 14 Feb 2020 15:49:11 +0200 Subject: [PATCH 13/13] modify namespace --- R-proj/NAMESPACE | 1 - 1 file changed, 1 deletion(-) diff --git a/R-proj/NAMESPACE b/R-proj/NAMESPACE index 36f46c881..317a83a2c 100644 --- a/R-proj/NAMESPACE +++ b/R-proj/NAMESPACE @@ -22,7 +22,6 @@ exportPattern("^[[:alpha:]]+") importFrom("methods","new") importFrom("stats", "cov") importFrom("utils","read.csv") -importFrom("stats", "cov") importFrom(Rcpp,evalCpp) importFrom(Rcpp,loadModule) useDynLib(volesti, .registration=TRUE)