Exact volume

R-proj/DESCRIPTION

@@ -15,7 +15,7 @@ Version: 1.1.0
 Date: 2020-01-23
 Maintainer: Vissarion Fisikopoulos <[email protected]>
 Depends: Rcpp (>= 0.12.17)
-Imports: methods
+Imports: methods, stats
 LinkingTo: Rcpp, RcppEigen, BH
 Suggests: testthat
 Encoding: UTF-8

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)

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.

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)

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)

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)
 

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)

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)
 }

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<Rcpp::Reference> P = R_NilValue, Rcpp::Nullable<std::string> body = R_NilValue,
-                 Rcpp::Nullable<Rcpp::List> Parameters = R_NilValue){
+double exact_vol(Rcpp::Nullable<Rcpp::Reference> P) {
 
     typedef double NT;
-    typedef Cartesian<NT>    Kernel;
-    typedef typename Kernel::Point    Point;
-    typedef boost::mt19937    RNGType;
-    typedef Eigen::Matrix<NT,Eigen::Dynamic,1> VT;
-    typedef Eigen::Matrix<NT,Eigen::Dynamic,Eigen::Dynamic> MT;
-
-    typedef Zonotope<Point> zonotope;
-    typedef VPolytope <Point, RNGType> Vpolytope;
-    zonotope ZP;
-    Vpolytope VP;
-    int type, dim;
-    NT vol, rad;
-
-    if (P.isNotNull()) {
-        type = Rcpp::as<Rcpp::Reference>(P).field("type");
-        if (!body.isNotNull()) {
-            dim = Rcpp::as<Rcpp::Reference>(P).field("dimension");
-            if (type == 3) {
-                ZP.init(dim, Rcpp::as<MT>(Rcpp::as<Rcpp::Reference>(P).field("G")),
-                        VT::Ones(Rcpp::as<MT>(Rcpp::as<Rcpp::Reference>(P).field("G")).rows()));
-                vol = exact_zonotope_vol<NT>(ZP);
-            } else if (type == 2) {
-                if (Rcpp::as<MT>(Rcpp::as<Rcpp::Reference>(P).field("V")).rows() == Rcpp::as<MT>(Rcpp::as<Rcpp::Reference>(P).field("V")).cols()+1) {
-                    MT V = Rcpp::as<MT>(Rcpp::as<Rcpp::Reference>(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<Rcpp::List>(Parameters).containsElementNamed("dimension")) {
-                    dim = Rcpp::as<int>(Rcpp::as<Rcpp::List>(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<std::string>(body).compare(std::string("simplex"))==0) {
-                vol = 1.0 / factorial(NT(dim));
-            } else if (Rcpp::as<std::string>(body).compare(std::string("cross"))==0) {
-                vol = std::pow(2.0, NT(dim));
-                vol = vol / factorial(NT(dim));
-            } else if (Rcpp::as<std::string>(body).compare(std::string("cube"))==0) {
-                vol = std::pow(2.0, NT(dim));
-            } else if (Rcpp::as<std::string>(body).compare(std::string("hypersphere"))==0) {
-                if (Rcpp::as<Rcpp::List>(Parameters).containsElementNamed("radius")) {
-                    rad = Rcpp::as<NT>(Rcpp::as<Rcpp::List>(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 <NT> Kernel;
+    typedef typename Kernel::Point Point;
+    typedef boost::mt19937 RNGType;
+    typedef Eigen::Matrix<NT, Eigen::Dynamic, 1> VT;
+    typedef Eigen::Matrix <NT, Eigen::Dynamic, Eigen::Dynamic> MT;
+
+    if (NT(Rcpp::as<Rcpp::Reference>(P).field("volume")) > 0.0) {
+        return NT(Rcpp::as<Rcpp::Reference>(P).field("volume"));
+    }
+
+    int type = Rcpp::as<Rcpp::Reference>(P).field("type"), dim;
+    NT vol;
+
+    if (type == 2) {
+
+        typedef VPolytope <Point, RNGType> Vpolytope;
+        Vpolytope VP;
+        dim = Rcpp::as<Rcpp::Reference>(P).field("dimension");
+
+        if (Rcpp::as<MT>(Rcpp::as<Rcpp::Reference>(P).field("V")).rows() ==
+            Rcpp::as<MT>(Rcpp::as<Rcpp::Reference>(P).field("V")).cols() + 1) {
+
+            MT V = Rcpp::as<MT>(Rcpp::as<Rcpp::Reference>(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 <Point> zonotope;
+        zonotope ZP;
+        dim = Rcpp::as<Rcpp::Reference>(P).field("dimension");
+
+        ZP.init(dim, Rcpp::as<MT>(Rcpp::as<Rcpp::Reference>(P).field("G")),
+                VT::Ones(Rcpp::as<MT>(Rcpp::as<Rcpp::Reference>(P).field("G")).rows()));
+        vol = exact_zonotope_vol<NT>(ZP);
+
+    } else {
+        throw Rcpp::exception("Volume unknown!");
     }
 
     return vol;