Marginal diagnostics and optimizations

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,7 +157,7 @@ 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.
 #'
@@ -168,8 +170,8 @@ geweke <- function(samples, frac_first = NULL, frac_last = NULL) {
 #' 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/gen_birkhoff.R

@@ -0,0 +1,27 @@
+#' Generator function for Birkhoff polytope
+#' 
+#' This function can be used to generate the full dimensional \eqn{n}-Birkhoff polytope in H-representation.
+#' 
+#' @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/src/geweke.cpp

@@ -27,6 +27,8 @@
 //' \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
 // [[Rcpp::export]]
 bool geweke(Rcpp::NumericMatrix samples, 
             Rcpp::Nullable<double> frac_first = R_NilValue, 

R-proj/src/inner_ball.cpp

@@ -20,7 +20,7 @@
 //' 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.
 //'
@@ -35,7 +35,7 @@
 //' @export
 // [[Rcpp::export]]
 Rcpp::NumericVector inner_ball(Rcpp::Reference P, 
-                               Rcpp::Nullable<std::string> method = R_NilValue) {
+                               Rcpp::Nullable<bool> lpsolve = R_NilValue) {
 
     typedef double NT;
     typedef Cartesian<NT>    Kernel;
@@ -50,33 +50,27 @@ Rcpp::NumericVector inner_ball(Rcpp::Reference P,
     unsigned int n = P.field("dimension"), type = P.field("type");
 
     std::pair <Point, NT> InnerBall;
-    std::string method_rcpp = (!method.isNotNull()) ? std::string("lpsolve") : Rcpp::as<std::string>(method);
+    bool lp_solve = (lpsolve.isNotNull()) ? Rcpp::as<bool>(lpsolve) : false;
 
     switch (type) {
         case 1: {
             // Hpolytope
             Hpolytope HP;
             HP.init(n, Rcpp::as<MT>(P.field("A")), Rcpp::as<VT>(P.field("b")));
-            if(method_rcpp.compare(std::string("lpsolve")) == 0) {
+            if (lp_solve) {
+                InnerBall = ComputeChebychevBall<NT, Point>(HP.get_mat(), HP.get_vec());
+            } else {
                 InnerBall = HP.ComputeInnerBall();
-                break;
-            } else if(method_rcpp.compare(std::string("ipm")) == 0) {
-                HP.normalize();
-                NT tol = 0.00000001;
-                std::pair<VT, NT> res = max_inscribed_ball(HP.get_mat(), HP.get_vec(), 150, tol);
-                InnerBall.second = res.second;
-                InnerBall.first = Point(res.first);
-            } else {         
-                throw Rcpp::exception("Unknown method!");
             }
-
+            if (InnerBall.second < 0.0) throw Rcpp::exception("Unable to compute a feasible point.");
             break;
         }
         case 2: {
             // Vpolytope
             Vpolytope VP;
             VP.init(n, Rcpp::as<MT>(P.field("V")), VT::Ones(Rcpp::as<MT>(P.field("V")).rows()));
             InnerBall = VP.ComputeInnerBall();
+            if (InnerBall.second < 0.0) throw Rcpp::exception("Unable to compute a feasible point.");
             break;
         }
         case 3: {
@@ -85,6 +79,7 @@ Rcpp::NumericVector inner_ball(Rcpp::Reference P,
             InnerBall = ZP.ComputeInnerBall();
             ZP.init(n, Rcpp::as<MT>(P.field("G")), VT::Ones(Rcpp::as<MT>(P.field("G")).rows()));
             InnerBall = ZP.ComputeInnerBall();
+            if (InnerBall.second < 0.0) throw Rcpp::exception("Unable to compute a feasible point.");
             break;
         }
         case 4: {
@@ -97,6 +92,7 @@ Rcpp::NumericVector inner_ball(Rcpp::Reference P,
             VPcVP.init(VP1, VP2);
             if (!VPcVP.is_feasible()) throw Rcpp::exception("Empty set!");
             InnerBall = VPcVP.ComputeInnerBall();
+            if (InnerBall.second < 0.0) throw Rcpp::exception("Unable to compute a feasible point.");
             break;
         }
     }

R-proj/src/poly_gen.cpp

@@ -101,6 +101,9 @@ Rcpp::NumericMatrix poly_gen (int kind_gen, bool Vpoly_gen, bool Zono_gen, int d
 
             case 6:
                 return extractMatPoly(random_hpoly<Hpolytope, RNGType>(dim_gen, m_gen, seed_rcpp));
+
+            case 7:
+                return extractMatPoly(gen_birk<Hpolytope>(dim_gen));
 
         }
     }

R-proj/src/psrf.cpp → R-proj/src/psrf_multivariate.cpp

@@ -14,7 +14,7 @@
 #include <boost/random/uniform_int.hpp>
 #include <boost/random/normal_distribution.hpp>
 #include <boost/random/uniform_real_distribution.hpp>
-#include "diagnostics/psrf.hpp"
+#include "diagnostics/multivariate_psrf.hpp"
 
 //' Gelman-Rubin Potential Scale Reduction Factor (PSRF)
 //'
@@ -27,12 +27,14 @@
 //' \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.
+//'
+//' @export
 // [[Rcpp::export]]
-double psrf(Rcpp::NumericMatrix samples)
+double psrf_multivariate(Rcpp::NumericMatrix samples)
 {
     typedef double NT;
     typedef Eigen::Matrix<NT,Eigen::Dynamic,1> VT;
     typedef Eigen::Matrix<NT,Eigen::Dynamic,Eigen::Dynamic> MT;    
 
-    return perform_psrf<NT, VT>(Rcpp::as<MT>(samples));
+    return multivariate_psrf<NT, VT>(Rcpp::as<MT>(samples));
 }

R-proj/src/psrf_univariate.cpp

@@ -0,0 +1,58 @@
+// [[Rcpp::depends(BH)]]
+
+// VolEsti (volume computation and sampling library)
+
+// Copyright (c) 20014-2020 Vissarion Fisikopoulos
+// Copyright (c) 2018-2020 Apostolos Chalkis
+
+//Contributed and/or modified by Alexandros Manochis, as part of Google Summer of Code 2020 program.
+
+#include <Rcpp.h>
+#include <RcppEigen.h>
+#include <chrono>
+#include <boost/random.hpp>
+#include <boost/random/uniform_int.hpp>
+#include <boost/random/normal_distribution.hpp>
+#include <boost/random/uniform_real_distribution.hpp>
+#include "diagnostics/univariate_psrf.hpp"
+#include "diagnostics/interval_psrf.hpp"
+
+//' 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
+// [[Rcpp::export]]
+Rcpp::NumericVector psrf_univariate(Rcpp::NumericMatrix samples,
+                                  Rcpp::Nullable<std::string> method = R_NilValue)
+{
+    typedef double NT;
+    typedef Eigen::Matrix<NT,Eigen::Dynamic,1> VT;
+    typedef Eigen::Matrix<NT,Eigen::Dynamic,Eigen::Dynamic> MT;
+
+    VT scores(samples.nrow());
+
+    std::string method_rcpp = std::string("normal");
+    if(method.isNotNull()) {
+        method_rcpp =  Rcpp::as<std::string>(method);
+    }
+
+    if (method_rcpp.compare(std::string("normal")) == 0) {
+        scores = univariate_psrf<NT, VT>(Rcpp::as<MT>(samples));
+    } else if(method_rcpp.compare(std::string("interval")) == 0) {
+        scores = interval_psrf<VT, NT>(Rcpp::as<MT>(samples));
+    } else {
+        throw Rcpp::exception("Unknown method!");
+    }
+
+    return Rcpp::wrap(scores);
+}

R-proj/src/raftery.cpp

@@ -28,8 +28,9 @@
 //' \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
 // [[Rcpp::export]]
-
 Rcpp::NumericMatrix raftery(Rcpp::NumericMatrix samples,
                             Rcpp::Nullable<double> q = R_NilValue,
                             Rcpp::Nullable<double> r = R_NilValue,