Marginal diagnostics and optimizations (#108)

* improve implementations

* improve computation of inner ball for H-polytopes

* improve inner ball computation and improve tests

* update termination criterions in rounding methods

* update Rd file of Birkhoff R generator

* improve univariate psrf implemetations

* fix c++ tests

* update parameters in rounding and minor improvements

* fix gcc tests and improve birkhoff generator Rd file

* fix c++ tests

* improve R tests, remove ine and ext files

* fix c++ tests

* fix c++ tests for clang

* fix c++ tests

* merge and improve R examples

* change r tests and examples for inner_ball function

* improve R examples

* change priority in inner ball computation

* use only lpsolve for inner ball computation

R-proj/NAMESPACE

@@ -6,6 +6,7 @@ export(direct_sampling)
 export(exact_vol)
 export(file_to_polytope)
 export(frustum_of_simplex)
+export(gen_birkhoff)
 export(gen_cross)
 export(gen_cube)
 export(gen_prod_simplex)
@@ -15,8 +16,12 @@ export(gen_rand_zonotope)
 export(gen_simplex)
 export(gen_skinny_cube)
 export(get_full_dimensional_polytope)
+export(geweke)
 export(inner_ball)
 export(loadSdpaFormatFile)
+export(psrf_multivariate)
+export(psrf_univariate)
+export(raftery)
 export(readSdpaFormatFile)
 export(rotate_polytope)
 export(round_polytope)

R-proj/R/RcppExports.R

@@ -145,6 +145,8 @@ full_dimensional_polytope <- function(P) {
 #' \dQuote{Evaluating the accuracy of sampling-based approaches to the calculation of posterior moments,} \emph{ In Bayesian Statistics 4. Proceedings of the Fourth Valencia International Meeting,} 1992.}
 #'
 #' @return A boolean to denote if the result of Geweke diagnostic: (i)  false if the null hypothesis is rejected, (ii) true if the null hypothesis is not rejected.
+#'
+#' @export
 geweke <- function(samples, frac_first = NULL, frac_last = NULL) {
     .Call(`_volesti_geweke`, samples, frac_first, frac_last)
 }
@@ -155,21 +157,21 @@ geweke <- function(samples, frac_first = NULL, frac_last = NULL) {
 #' For both zonotopes and V-polytopes the function computes the minimum \eqn{r} s.t.: \eqn{ r e_i \in P} for all \eqn{i=1, \dots ,d}. Then the ball centered at the origin with radius \eqn{r/ \sqrt{d}} is an inscribed ball.
 #'
 #' @param P A convex polytope. It is an object from class (a) Hpolytope or (b) Vpolytope or (c) Zonotope or (d) VpolytopeIntersection.
-#' @param method Optional. A string to declare the method to be used: (i) \code{'lpsolve'} to use lpsolve library, (ii) \code{'ipm'} to use an interior point method which solves the corresponding linear program. The default method is \code{'lpsolve'}.
+#' @param lpsolve Optional. A boolean variable to compute the Chebychev ball of an H-polytope using the lpsolve library.
 #'
 #' @return A \eqn{(d+1)}-dimensional vector that describes the inscribed ball. The first \eqn{d} coordinates corresponds to the center of the ball and the last one to the radius.
 #'
 #' @examples
 #' # compute the Chebychev ball of the 2d unit simplex
-#' P = gen_simplex(2,'H')
+#' P = gen_cube(10,'H')
 #' ball_vec = inner_ball(P)
 #'
 #' # compute an inscribed ball of the 3-dimensional unit cube in V-representation
 #' P = gen_cube(3, 'V')
 #' ball_vec = inner_ball(P)
 #' @export
-inner_ball <- function(P, method = NULL) {
-    .Call(`_volesti_inner_ball`, P, method)
+inner_ball <- function(P, lpsolve = NULL) {
+    .Call(`_volesti_inner_ball`, P, lpsolve)
 }
 
 #' An internal Rccp function as a polytope generator
@@ -199,8 +201,28 @@ poly_gen <- function(kind_gen, Vpoly_gen, Zono_gen, dim_gen, m_gen, seed = NULL)
 #' \dQuote{General Methods for Monitoring Convergence of Iterative Simulations,} \emph{Journal of Computational and Graphical Statistics,} 1998.}
 #'
 #' @return The value of multivariate PSRF by S. Brooks and A. Gelman.
-psrf <- function(samples) {
-    .Call(`_volesti_psrf`, samples)
+#'
+#' @export
+psrf_multivariate <- function(samples) {
+    .Call(`_volesti_psrf_multivariate`, samples)
+}
+
+#' Gelman-Rubin and Brooks-Gelman Potential Scale Reduction Factor (PSRF) for each marginal
+#'
+#' @param samples A matrix that contans column-wise the sampled points from a geometric random walk.
+#' @param method A string to reauest diagnostic: (i) \code{'normal'} for psrf of Gelman-Rubin and (ii) \code{'interval'} for psrf of Brooks-Gelman.
+#'
+#' @references \cite{Gelman, A. and Rubin, D. B.,
+#' \dQuote{Inference from iterative simulation using multiple sequences,} \emph{Statistical Science,} 1992.}
+#'
+#' @references \cite{Brooks, S. and Gelman, A.,
+#' \dQuote{General Methods for Monitoring Convergence of Iterative Simulations,} \emph{Journal of Computational and Graphical Statistics,} 1998.}
+#'
+#' @return A vector that contains the values of PSRF for each coordinate
+#'
+#' @export
+psrf_univariate <- function(samples, method = NULL) {
+    .Call(`_volesti_psrf_univariate`, samples, method)
 }
 
 #' Raftery and Lewis MCMC diagnostic
@@ -214,6 +236,8 @@ psrf <- function(samples) {
 #' \dQuote{How many iterations in the Gibbs sampler?,} \emph{Bayesian Statistics 4. Proceedings of the Fourth Valencia International Meeting,} 1992.}
 #'
 #' @return (i) The number of draws required for burn-in, (ii) the skip parameter for 1st-order Markov chain, (iii) the skip parameter sufficient to get independence chain, (iv) the number of draws required to achieve r precision, (v) the number of draws if the chain is white noise, (vi) the I-statistic from Raftery and Lewis (1992).
+#'
+#' @export
 raftery <- function(samples, q = NULL, r = NULL, s = NULL) {
     .Call(`_volesti_raftery`, samples, q, r, s)
 }

R-proj/R/file_to_polytope.R

@@ -8,10 +8,6 @@
 #'
 #' @return A polytope class. If the path corresponds to an ine file then the return value represents a H-polytope. If it corresponds to an ext file the return value represents a V-polytope (default choice) or a zonotope if the second argument is TRUE.
 #'
-#' @examples
-#' # give the path to birk4.ine
-#' path = system.file('extdata', package = 'volesti')
-#' P = file_to_polytope(paste0(path,'/birk4.ine'))
 #' @export
 #' @useDynLib volesti, .registration=TRUE
 #' @importFrom Rcpp evalCpp
@@ -20,7 +16,7 @@
 #' @importFrom "stats" "cov"
 #' @importFrom "methods" "new"
 #' @exportPattern "^[[:alpha:]]+"
-file_to_polytope <- function(path, zonotope){
+file_to_polytope <- function(path, zonotope = FALSE){
 
   ineorext=substr(path, start = nchar(path) - 2, stop = nchar(path))
   if(ineorext!="ine" && ineorext!="ext") {

R-proj/R/gen_birkhoff.R

@@ -0,0 +1,28 @@
+#' Generator function for Birkhoff polytope
+#' 
+#' This function can be used to generate the full dimensional \eqn{n}-Birkhoff polytope in H-representation.
+#' The dimension of the generated polytope is \eqn{(n-1)^2}.
+#' 
+#' @param n The order of the Birkhoff polytope
+#' 
+#' @return A polytope class representing the full dimensional \eqn{n}-Birkhoff polytope in H-representation.
+#' @examples 
+#' # generate the Birkhoff polytope of order 5
+#' P = gen_birkhoff(5)
+#' @export
+gen_birkhoff <- function(n) {
+
+  kind_gen = 7
+  m_gen = 0
+
+  Mat = poly_gen(kind_gen, FALSE, FALSE, n, m_gen)
+
+  # first column is the vector b
+  b = Mat[,1]
+  Mat = Mat[,-c(1)]
+
+  P = Hpolytope$new(Mat, b)
+
+  return(P)
+
+}

R-proj/R/gen_cross.R

@@ -38,4 +38,4 @@ gen_cross <- function(dimension, representation) {
   }
 
   return(P)
-}
+}

R-proj/R/round_polytope.R

@@ -19,17 +19,15 @@
 #' 
 #'
 #' @examples
-#' # rotate a H-polytope (2d unit simplex)
-#' 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)
+#' # round a 5d skinny cube
+#' P = gen_skinny_cube(5)
 #' listHpoly = round_polytope(P)
 #' 
-#' # rotate a V-polytope (3d unit cube) using Random Directions HnR with step equal to 50
+#' # round a V-polytope (3d unit cube)
 #' P = gen_cube(3, 'V')
 #' ListVpoly = round_polytope(P)
 #' 
-#' # round a 2-dimensional zonotope defined by 6 generators using ball walk
+#' # round a 2-dimensional zonotope defined by 6 generators
 #' Z = gen_rand_zonotope(2,6)
 #' ListZono = round_polytope(Z)
 #' @export

R-proj/R/zonotope_approximation.R

@@ -20,7 +20,7 @@
 #' 
 #' @examples
 #' # over-approximate a 2-dimensional zonotope with 10 generators and compute the ratio of fitness
-#' Z = gen_rand_zonotope(2,8)
+#' Z = gen_rand_zonotope(2,12)
 #' retList = zonotope_approximation(Z = Z)
 #' 
 #' @export
@@ -40,4 +40,3 @@ zonotope_approximation <- function(Z, fit_ratio = NULL, settings = NULL, seed =
   return(PP)
 
 }
-