From 3ef59d043afd7f490ddfb86e033f8f1b5d6d51a5 Mon Sep 17 00:00:00 2001 From: Tolis Date: Sun, 16 Feb 2020 20:36:28 +0200 Subject: [PATCH 01/14] use cdhr in rounding, improve t-test iterations, change diameter of H-polytopes, minor mpdifications in both interfaces. --- R-proj/R/round_polytope.R | 4 +--- include/annealing/ball_annealing.h | 9 +++++---- include/convex_bodies/hpolytope.h | 4 ++-- include/rounding.h | 7 +++++-- include/volume/volume.h | 6 +++++- test/vol.cpp | 18 +++++++++++++----- 6 files changed, 31 insertions(+), 17 deletions(-) diff --git a/R-proj/R/round_polytope.R b/R-proj/R/round_polytope.R index f980b8639..16b8367ae 100644 --- a/R-proj/R/round_polytope.R +++ b/R-proj/R/round_polytope.R @@ -45,10 +45,8 @@ round_polytope <- function(P, random_walk = NULL, walk_length = NULL, parameters PP = list("P" = Vpolytope$new(A), "round_value" = ret_list$round_value) }else if (type == 3) { PP = list("P" = Zonotope$new(A), "round_value" = ret_list$round_value) - } else if(type == 1) { - PP = list("P" = Hpolytope$new(A,b), "round_value" = ret_list$round_value) } else { - PP = list("P" = VPolyintersectVPoly$new(V1 = t(Mat[,dim(P$V1)[1]]), V2 = t(Mat[,dim(P$V2)[1]])), "round_value" = ret_list$round_value) + PP = list("P" = Hpolytope$new(A,b), "round_value" = ret_list$round_value) } return(PP) } diff --git a/include/annealing/ball_annealing.h b/include/annealing/ball_annealing.h index e8cbf82eb..33ad1cc5f 100644 --- a/include/annealing/ball_annealing.h +++ b/include/annealing/ball_annealing.h @@ -16,14 +16,14 @@ bool check_convergence(ConvexBody &P, PointList &randPoints, const NT &lb, const std::vector ratios; std::pair mv; - int m = randPoints.size()/nu; + int m = randPoints.size()/nu, i = 1; NT T, rs, alpha_check = 0.01; size_t countsIn = 0; - for(typename PointList::iterator pit=randPoints.begin(); pit!=randPoints.end(); ++pit){ + for(typename std::list::iterator pit=randPoints.begin(); pit!=randPoints.end(); ++pit, i++){ if (P.is_in(*pit)==-1) countsIn++; - if ((std::distance(randPoints.begin(), pit)+1) % m == 0) { + if (i % m == 0) { ratios.push_back(NT(countsIn)/m); countsIn = 0; if (ratios.size()>1 && precheck) { @@ -47,7 +47,8 @@ bool check_convergence(ConvexBody &P, PointList &randPoints, const NT &lb, const ratio = mv.first; rs = std::sqrt(mv.second); boost::math::students_t dist(nu - 1); - T = rs*(boost::math::quantile(boost::math::complement(dist, alpha)) / std::sqrt(NT(nu))); + T = rs*(boost::math::quantile(boost::math::complement(dist, alpha)) + / std::sqrt(NT(nu))); if (ratio > lb + T) { if (lastball) return true; if ((precheck && ratio < ub - T) || (!precheck && ratio < ub + T)) return true; diff --git a/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h index 12ae42bff..487358188 100644 --- a/include/convex_bodies/hpolytope.h +++ b/include/convex_bodies/hpolytope.h @@ -140,9 +140,9 @@ class HPolytope{ void comp_diam(NT &diam, const NT &cheb_rad) { if(cheb_rad < 0.0) { - diam = 2.0 * std::sqrt(NT(_d)) * ComputeInnerBall().second; + diam = 4.0 * std::sqrt(NT(_d)) * ComputeInnerBall().second; } else { - diam = 2.0 * std::sqrt(NT(_d)) * cheb_rad; + diam = 4.0 * std::sqrt(NT(_d)) * cheb_rad; } } diff --git a/include/rounding.h b/include/rounding.h index d79fd82a4..65d6a6389 100644 --- a/include/rounding.h +++ b/include/rounding.h @@ -34,11 +34,14 @@ std::pair rounding_min_ellipsoid(Polytope &P , const std::pair::iterator pit=randPoints.begin(); pit!=randPoints.end(); ++pit){ current_dist=(*pit-c).squared_length(); diff --git a/include/volume/volume.h b/include/volume/volume.h index 906f711e2..8ba24da4d 100644 --- a/include/volume/volume.h +++ b/include/volume/volume.h @@ -96,7 +96,7 @@ NT volume(Polytope &P, Point p = get_point_on_Dsphere(n, radius); std::list randPoints; //ds for storing rand points //use a large walk length e.g. 1000 - + rand_point_generator(P, p, 1, 50*n, randPoints, var); double tstart2 = (double)clock()/(double)CLOCKS_PER_SEC; @@ -207,6 +207,10 @@ NT volume(Polytope &P, < Date: Mon, 17 Feb 2020 10:45:03 +0200 Subject: [PATCH 02/14] improve c++ interface, add c++ test for cooling bodies with billiard walk --- include/volume/volume.h | 4 - test/CMakeLists.txt | 9 + test/cooling_bodies_bill_test.cpp | 282 ++++++++++++++++++++++++++++++ test/vol.cpp | 44 ++--- 4 files changed, 314 insertions(+), 25 deletions(-) create mode 100644 test/cooling_bodies_bill_test.cpp diff --git a/include/volume/volume.h b/include/volume/volume.h index 8ba24da4d..f5aeac3a3 100644 --- a/include/volume/volume.h +++ b/include/volume/volume.h @@ -207,10 +207,6 @@ NT volume(Polytope &P, <) add_executable (VpolyVol_test VpolyVol_test.cpp $) add_executable (ZonotopeVol_test ZonotopeVol_test.cpp $) + add_executable (cool_bodies_bill_test cooling_bodies_bill_test.cpp $) #add_executable (ZonotopeVolCV_test ZonotopeVolCV_test.cpp $) add_test(NAME volume_cube COMMAND volume_test -tc=cube) @@ -108,6 +109,13 @@ else () add_test(NAME cheb_simplex COMMAND cheb_test -tc=cheb_simplex) add_test(NAME cheb_skinny_cube COMMAND cheb_test -tc=cheb_skinny_cube) + add_test(NAME cool_bodies_cube COMMAND cool_bodies_bill_test -tc=cube) + add_test(NAME cool_bodies_cross COMMAND cool_bodies_bill_test -tc=cross) + add_test(NAME cool_bodies_birkhoff COMMAND cool_bodies_bill_test -tc=birk) + add_test(NAME cool_bodies_prod_simplex COMMAND cool_bodies_bill_test -tc=prod_simplex) + add_test(NAME cool_bodies_simplex COMMAND cool_bodies_bill_test -tc=simplex) + add_test(NAME cool_bodies_skinny_cube COMMAND cool_bodies_bill_test -tc=skinny_cube) + #add_test(NAME round_skinny_cube COMMAND rounding_test -tc=round_skinny_cube) #add_test(NAME round_rot_skinny_cube COMMAND rounding_test -tc=round_rot_skinny_cube) @@ -121,6 +129,7 @@ else () TARGET_LINK_LIBRARIES(VpolyCV_test ${LP_SOLVE}) TARGET_LINK_LIBRARIES(VpolyVol_test ${LP_SOLVE}) TARGET_LINK_LIBRARIES(ZonotopeVol_test ${LP_SOLVE}) + TARGET_LINK_LIBRARIES(cool_bodies_bill_test ${LP_SOLVE}) #TARGET_LINK_LIBRARIES(ZonotopeVolCV_test ${LP_SOLVE}) endif() diff --git a/test/cooling_bodies_bill_test.cpp b/test/cooling_bodies_bill_test.cpp new file mode 100644 index 000000000..4a9b2331e --- /dev/null +++ b/test/cooling_bodies_bill_test.cpp @@ -0,0 +1,282 @@ +#include "doctest.h" +#include +#include "Eigen/Eigen" +#define VOLESTI_DEBUG +#include +#include "random.hpp" +#include "random/uniform_int.hpp" +#include "random/normal_distribution.hpp" +#include "random/uniform_real_distribution.hpp" +#include +#include +#include +#include +#include +#include +#include "cartesian_geom/cartesian_kernel.h" +#include "vars.h" +#include "hpolytope.h" +#include "vpolytope.h" +#include "zpolytope.h" +#include "ball.h" +#include "ballintersectconvex.h" +#include "vpolyintersectvpoly.h" +#include "samplers.h" +#include "rounding.h" +#include "rotating.h" +#include "linear_extensions.h" +#include "gaussian_annealing.h" +#include "cooling_balls.h" +#include "cooling_hpoly.h" +#include "misc.h" +#include "polytope_generators.h" +#include + +template +NT factorial(NT n) +{ + return (n == 1 || n == 0) ? 1 : factorial(n - 1) * n; +} + +template +void test_cool_bodies(Polytope &HP, NT expected, NT tolerance=0.1, bool round = false, NT diam = -1.0) +{ + + typedef typename Polytope::PolytopePoint Point; + + // Setup the parameters + int n = HP.dimension(); + int walk_len=10 + n/10; + int nexp=1, n_threads=1; + NT e=0.1, err=0.0000000001, diameter = diam, round_val = 1.0; + int rnum = std::pow(e,-2) * 400 * n * std::log(n); + unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); + RNGType rng(seed); + boost::normal_distribution<> rdist(0,1); + boost::random::uniform_real_distribution<>(urdist); + boost::random::uniform_real_distribution<> urdist1(-1,1); + + + std::pair InnerBall; + if (round) { + InnerBall = HP.ComputeInnerBall(); + vars var2(rnum,n,walk_len,n_threads,err,e,0,0,0,0,0.0,rng, + urdist,urdist1,-1.0,false,false,false,false,false,false,false,false,true); + std::pair res_round = rounding_min_ellipsoid(HP , InnerBall, var2); + round_val = res_round.first; + } + + InnerBall = HP.ComputeInnerBall(); + if(diameter < 0.0) { + HP.comp_diam(diameter, InnerBall.second); + } + + vars var(rnum,n,walk_len,n_threads,err,e,0,0,0,InnerBall.second,diameter,rng, + urdist,urdist1,-1.0,false,false,false,false,false,false,false,false,true); + + NT lb = 0.1, ub = 0.15, p = 0.75, rmax = 0.0, alpha = 0.2; + int W = 250, NNu = 140, nu =10; + bool win2 = false; + vars_ban var_ban(lb, ub, p, rmax, alpha, W, NNu, nu, win2); + + //Compute chebychev ball// + std::pair CheBall; + + // Estimate the volume + std::cout << "Number type: " << typeid(NT).name() << std::endl; + NT vol = 0; + unsigned int const num_of_exp = 10; + for (unsigned int i=0; i(10, false); + test_cool_bodies(P, 1024.0); + + std::cout << "--- Testing volume of H-cube20" << std::endl; + P = gen_cube(20, false); + test_cool_bodies(P, 1048576.0); + + std::cout << "--- Testing volume of H-cube30" << std::endl; + P = gen_cube(30, false); + test_cool_bodies(P, 1073742000.0); +} + +template +void call_test_cross(){ + typedef Cartesian Kernel; + typedef typename Kernel::Point Point; + typedef boost::mt19937 RNGType; + typedef HPolytope Hpolytope; + + std::cout << "--- Testing volume of H-cross10" << std::endl; + Hpolytope P = gen_cross(10, false); + test_cool_bodies(P, 0.0002821869); +} + +template +void call_test_birk() { + typedef Cartesian Kernel; + typedef typename Kernel::Point Point; + typedef boost::mt19937 RNGType; + typedef HPolytope Hpolytope; + Hpolytope P; + + std::cout << "--- Testing volume of H-birk3" << std::endl; + std::ifstream inp; + std::vector > Pin; + inp.open("../R-proj/inst/extdata/birk3.ine",std::ifstream::in); + read_pointset(inp,Pin); + P.init(Pin); + test_cool_bodies(P, 0.125); + + std::cout << "--- Testing volume of H-birk4" << std::endl; + std::ifstream inp2; + std::vector > Pin2; + inp2.open("../R-proj/inst/extdata/birk4.ine",std::ifstream::in); + read_pointset(inp2,Pin2); + P.init(Pin2); + test_cool_bodies(P, 0.000970018, 0.2); + + std::cout << "--- Testing volume of H-birk5" << std::endl; + std::ifstream inp3; + std::vector > Pin3; + inp3.open("../R-proj/inst/extdata/birk5.ine",std::ifstream::in); + read_pointset(inp3,Pin3); + P.init(Pin3); + test_cool_bodies(P, 0.000000225, 0.2); + + //std::cout << "--- Testing volume of H-birk6" << std::endl; + //std::ifstream inp4; + //std::vector > Pin4; + //inp4.open("../R-proj/inst/extdata/birk6.ine",std::ifstream::in); + //read_pointset(inp4,Pin4); + //P.init(Pin4); + //test_volume(P, 0.0000000000009455459196, 0.5); +} + +template +void call_test_prod_simplex() { + typedef Cartesian Kernel; + typedef typename Kernel::Point Point; + typedef boost::mt19937 RNGType; + typedef HPolytope Hpolytope; + Hpolytope P; + + std::cout << "--- Testing volume of H-prod_simplex5" << std::endl; + P = gen_prod_simplex(5); + test_cool_bodies(P, std::pow(1.0 / factorial(5.0), 2.0),0.15, true); + + std::cout << "--- Testing volume of H-prod_simplex10" << std::endl; + P = gen_prod_simplex(10); + test_cool_bodies(P, std::pow(1.0 / factorial(10.0), 2.0), 0.15, true); + + //std::cout << "--- Testing volume of H-prod_simplex15" << std::endl; + //P = gen_prod_simplex(15); + //test_volume(P, std::pow(1.0 / factorial(15.0), 2)); + + std::cout << "--- Testing volume of H-prod_simplex20" << std::endl; + P = gen_prod_simplex(20); + test_cool_bodies(P, std::pow(1.0 / factorial(20.0), 2.0), 0.25, true); +} + +template +void call_test_simplex() { + typedef Cartesian Kernel; + typedef typename Kernel::Point Point; + typedef boost::mt19937 RNGType; + typedef HPolytope Hpolytope; + Hpolytope P; + + std::cout << "--- Testing volume of H-simplex10" << std::endl; + P = gen_simplex(10, false); + test_cool_bodies(P, 1.0 / factorial(10.0), 0.1, false, 1.5); + + std::cout << "--- Testing volume of H-simplex20" << std::endl; + P = gen_simplex(20, false); + test_cool_bodies(P, 1.0 / factorial(20.0), 0.1, false, 1.5); + + std::cout << "--- Testing volume of H-simplex30" << std::endl; + P = gen_simplex(30, false); + test_cool_bodies(P, 1.0 / factorial(30.0), 0.1, false, 1.5); + + //std::cout << "--- Testing volume of H-simplex40" << std::endl; + //P = gen_simplex(40, false); + //test_volume(P, 1.0 / factorial(40.0)); + + //std::cout << "--- Testing volume of H-simplex50" << std::endl; + //P = gen_simplex(50, false); + //test_volume(P, 1.0 / factorial(50.0)); +} + +template +void call_test_skinny_cube() { + typedef Cartesian Kernel; + typedef typename Kernel::Point Point; + typedef boost::mt19937 RNGType; + typedef HPolytope Hpolytope; + Hpolytope P; + + std::cout << "--- Testing volume of H-skinny_cube10" << std::endl; + P = gen_skinny_cube(10); + test_cool_bodies(P, 102400.0, 0.1, true); + + //std::cout << "--- Testing volume of H-skinny_cube20" << std::endl; + //P = gen_skinny_cube(20); + //test_volume(P, 104857600.0); +} + + +TEST_CASE("cube") { +call_test_cube(); +//call_test_cube(); +//call_test_cube(); +} + +TEST_CASE("cross") { +call_test_cross(); +//call_test_cross(); +//call_test_cross(); +} + +TEST_CASE("birk") { +call_test_birk(); +//call_test_birk(); +//call_test_birk(); +} + +TEST_CASE("prod_simplex") { +call_test_prod_simplex(); +//call_test_prod_simplex(); +//call_test_prod_simplex(); +} + +TEST_CASE("simplex") { +call_test_simplex(); +//call_test_simplex(); +//call_test_simplex(); +} + +TEST_CASE("skinny_cube") { +call_test_skinny_cube(); +//call_test_skinny_cube(); +//call_test_skinny_cube(); +} + diff --git a/test/vol.cpp b/test/vol.cpp index 45520d9b1..b033fd606 100644 --- a/test/vol.cpp +++ b/test/vol.cpp @@ -418,12 +418,35 @@ int main(const int argc, const char** argv) return 0; } + if (!set_algo) { + if (Zono || Vpoly) { + CB = true; + } else { + if (n <= 200) { + CB = true; + } else { + CG = true; + } + } + } else { + if (!CB && !CG) { + if (!set_error) { + e = 1.0; + error = 1.0; + } + } + } + if (!user_randwalk) { if (Zono || Vpoly) { billiard = true; } else { cdhr = true; } + } else if (!CB && !CG && billiard) { + std::cout<<"Billiard is not supported for SOB algorithm. CDHR is used."< Date: Mon, 17 Feb 2020 14:02:35 +0200 Subject: [PATCH 03/14] improve new c++ test --- R-proj/src/volume.cpp | 8 +------- test/cooling_bodies_bill_test.cpp | 5 +++-- 2 files changed, 4 insertions(+), 9 deletions(-) diff --git a/R-proj/src/volume.cpp b/R-proj/src/volume.cpp index 0ef136cc5..b8fda75dc 100644 --- a/R-proj/src/volume.cpp +++ b/R-proj/src/volume.cpp @@ -203,6 +203,7 @@ double volume (Rcpp::Reference P, Rcpp::Nullable walk_length = R_ if (!random_walk.isNotNull()) { if ( type == 1 ){ cdhr = true; + win_len = 3*n*n+400; } else { billiard = true; win_len = 150; @@ -257,11 +258,6 @@ double volume (Rcpp::Reference P, Rcpp::Nullable walk_length = R_ CB = true; e = (!error.isNotNull()) ? 0.1 : Rcpp::as(error); walkL = (!walk_length.isNotNull()) ? 1 : Rcpp::as(walk_length); - win_len = (cdhr) ? 3*n*n+400 : 2*n*n+250; - if (billiard){ - win_len = 150; - NN = 125; - } } else { throw Rcpp::exception("Unknown method!"); @@ -319,8 +315,6 @@ double volume (Rcpp::Reference P, Rcpp::Nullable walk_length = R_ } } - //std::cout<<"Wlen = "<(20); - test_cool_bodies(P, std::pow(1.0 / factorial(20.0), 2.0), 0.25, true); + test_cool_bodies(P, std::pow(1.0 / factorial(20.0), 2.0), 0.2, true); } template From d8284a3a701347d95e6db35ce8049800d381c2b5 Mon Sep 17 00:00:00 2001 From: tolischal Date: Mon, 17 Feb 2020 15:36:34 +0200 Subject: [PATCH 04/14] add random generators --- include/convex_bodies/vpolytope.h | 2 + include/convex_bodies/zpolytope.h | 2 + include/generators/h_polytopes_gen.h | 47 +++++ ...nerators.h => known_polytope_generators.h} | 101 +--------- include/generators/v_polytopes_gen.h | 149 +++++++++++++++ include/generators/z_polytopes_gen.h | 123 ++++++++++++ test/VpolyCV_test.cpp | 2 +- test/VpolyVol_test.cpp | 2 +- test/ZonotopeVolCV_test.cpp | 4 +- test/ZonotopeVol_test.cpp | 4 +- test/chebychev_test.cpp | 2 +- test/cooling_bodies_bill_test.cpp | 2 +- test/generator.cpp | 176 +++++++++++++++--- test/vol.cpp | 19 +- test/volumeCV_test.cpp | 2 +- test/volume_test.cpp | 2 +- 16 files changed, 501 insertions(+), 138 deletions(-) create mode 100644 include/generators/h_polytopes_gen.h rename include/generators/{polytope_generators.h => known_polytope_generators.h} (76%) create mode 100644 include/generators/v_polytopes_gen.h create mode 100644 include/generators/z_polytopes_gen.h diff --git a/include/convex_bodies/vpolytope.h b/include/convex_bodies/vpolytope.h index e03ed1015..f88cb039d 100644 --- a/include/convex_bodies/vpolytope.h +++ b/include/convex_bodies/vpolytope.h @@ -141,6 +141,8 @@ class VPolytope{ conv_comb = (REAL *) malloc(Pin.size() * sizeof(*conv_comb)); conv_comb2 = (REAL *) malloc(Pin.size() * sizeof(*conv_comb2)); colno = (int *) malloc((V.rows()+1) * sizeof(*colno)); + colno_mem = (int *) malloc(V.rows() * sizeof(*colno_mem)); + conv_mem = (REAL *) malloc(V.rows() * sizeof(*conv_mem)); row = (REAL *) malloc((V.rows()+1) * sizeof(*row)); } diff --git a/include/convex_bodies/zpolytope.h b/include/convex_bodies/zpolytope.h index 98fef55d7..2f1021ab1 100644 --- a/include/convex_bodies/zpolytope.h +++ b/include/convex_bodies/zpolytope.h @@ -202,6 +202,8 @@ class Zonotope { conv_comb = (REAL *) malloc(Pin.size() * sizeof(*conv_comb)); colno = (int *) malloc((V.rows()+1) * sizeof(*colno)); row = (REAL *) malloc((V.rows()+1) * sizeof(*row)); + colno_mem = (int *) malloc((V.rows()) * sizeof(*colno_mem)); + row_mem = (REAL *) malloc((V.rows()) * sizeof(*row_mem)); compute_eigenvectors(V.transpose()); } diff --git a/include/generators/h_polytopes_gen.h b/include/generators/h_polytopes_gen.h new file mode 100644 index 000000000..047d77c83 --- /dev/null +++ b/include/generators/h_polytopes_gen.h @@ -0,0 +1,47 @@ +// VolEsti (volume computation and sampling library) + +// Copyright (c) 20012-2018 Vissarion Fisikopoulos +// Copyright (c) 2018 Apostolos Chalkis + +// Licensed under GNU LGPL.3, see LICENCE file + +#ifndef H_POLYTOPES_GEN_H +#define H_POLYTOPES_GEN_H + +#include +#include "samplers.h" + +template +Polytope random_hpoly(unsigned int dim, unsigned int m) { + + typedef typename Polytope::MT MT; + typedef typename Polytope::VT VT; + typedef typename Polytope::NT NT; + typedef typename Polytope::PolytopePoint Point; + + unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); + RNGType rng(seed); + boost::random::uniform_real_distribution<> urdist1(-10, 10); + Point p(dim); + typename std::vector::iterator pit; + MT A(m, dim); + VT b(m); + unsigned int j; + + for(unsigned int i=0; i(dim); + pit = p.iter_begin(); + j = 0; + for ( ; pit!=p.iter_end(); ++pit, ++j) { + A(i,j) = *pit; + } + b(i) = 10.0; + + } + Polytope HP; + HP.init(dim, A, b); + + return HP; +} + +#endif diff --git a/include/generators/polytope_generators.h b/include/generators/known_polytope_generators.h similarity index 76% rename from include/generators/polytope_generators.h rename to include/generators/known_polytope_generators.h index 9653f95c0..df52215ba 100644 --- a/include/generators/polytope_generators.h +++ b/include/generators/known_polytope_generators.h @@ -5,8 +5,8 @@ // Licensed under GNU LGPL.3, see LICENCE file -#ifndef POLYTOPE_GENERATORS_H -#define POLYTOPE_GENERATORS_H +#ifndef KNOWN_POLYTOPE_GENERATORS_H +#define KNOWN_POLYTOPE_GENERATORS_H #include #include "samplers.h" @@ -306,103 +306,6 @@ Polytope gen_skinny_cube(const unsigned int &dim, bool Vpoly = false) { return P; } - -template -Polytope gen_zonotope(const unsigned int &dim, const unsigned int &m) { - - typedef typename Polytope::MT MT; - typedef typename Polytope::VT VT; - typedef typename Polytope::NT NT; - - unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); - RNGType rng(seed); - boost::normal_distribution<> rdist(0, 1); - - MT A; - VT b; - A.resize(m, dim); - b.resize(m); - Polytope P; - - for (unsigned int i = 0; i < m; ++i) { - b(i) = 1.0; - for (unsigned int j = 0; j < dim; ++j) { - A(i,j) = rdist(rng); - } - } - - P.init(dim, A, b); - return P; -} - -template -Polytope random_vpoly(const unsigned int &dim, const unsigned int &k) { - - typedef typename Polytope::MT MT; - typedef typename Polytope::VT VT; - typedef typename Polytope::NT NT; - typedef typename Polytope::PolytopePoint Point; - - unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); - RNGType rng(seed); - - Point p; - typename std::vector::iterator pit; - MT V(k, dim); - unsigned int j; - - for (unsigned int i = 0; i < k; ++i) { - p = get_direction(dim); - pit = p.iter_begin(); - j = 0; - for ( ; pit!=p.iter_end(); ++pit, ++j) { - V(i,j) = *pit; - } - } - - Polytope VP; - VT b = VT::Ones(k); - VP.init(dim, V, b); - - return VP; - -} - -template -Polytope random_hpoly(const unsigned int &dim, const unsigned int &m) { - - typedef typename Polytope::MT MT; - typedef typename Polytope::VT VT; - typedef typename Polytope::NT NT; - typedef typename Polytope::PolytopePoint Point; - - unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); - RNGType rng(seed); - boost::random::uniform_real_distribution<> urdist1(-10, 10); - Point p(dim); - typename std::vector::iterator pit; - MT A(m, dim); - VT b(m); - unsigned int j; - - for(unsigned int i=0; i(dim); - pit = p.iter_begin(); - j = 0; - for ( ; pit!=p.iter_end(); ++pit, ++j) { - A(i,j) = *pit; - } - b(i) = 10.0; - - } - Polytope HP; - HP.init(dim, A, b); - - return HP; -} - - - /* * ToDo: brkhoff polytope generator template diff --git a/include/generators/v_polytopes_gen.h b/include/generators/v_polytopes_gen.h new file mode 100644 index 000000000..73e148bbf --- /dev/null +++ b/include/generators/v_polytopes_gen.h @@ -0,0 +1,149 @@ +// VolEsti (volume computation and sampling library) + +// Copyright (c) 20012-2018 Vissarion Fisikopoulos +// Copyright (c) 2018 Apostolos Chalkis + +// Licensed under GNU LGPL.3, see LICENCE file + +#ifndef V_POLYTOPES_GEN_H +#define V_POLYTOPES_GEN_H + +#include +#include "samplers.h" + + +template +void removeRow(MT &matrix, unsigned int rowToRemove) +{ + unsigned int numRows = matrix.rows()-1; + unsigned int numCols = matrix.cols(); + + if( rowToRemove < numRows ) + matrix.block(rowToRemove,0,numRows-rowToRemove,numCols) = matrix.bottomRows(numRows-rowToRemove); + + matrix.conservativeResize(numRows,numCols); +} + +template +Polytope random_vpoly(unsigned int dim, unsigned int k) { + + typedef typename Polytope::MT MT; + typedef typename Polytope::VT VT; + typedef typename Polytope::NT NT; + typedef typename Polytope::PolytopePoint Point; + + unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); + RNGType rng(seed); + + Point p; + typename std::vector::iterator pit; + MT V(k, dim); + unsigned int j; + + for (unsigned int i = 0; i < k; ++i) { + p = get_direction(dim); + pit = p.iter_begin(); + j = 0; + for ( ; pit!=p.iter_end(); ++pit, ++j) { + V(i,j) = *pit; + } + } + + Polytope VP; + VT b = VT::Ones(k); + VP.init(dim, V, b); + + return VP; + +} + + +template +Polytope random_vpoly_incube(unsigned int d, unsigned int k) { + + typedef typename Polytope::MT MT; + typedef typename Polytope::VT VT; + typedef typename Polytope::NT NT; + typedef typename Polytope::PolytopePoint Point; + + REAL *conv_mem; + int *colno_mem; + + conv_mem = (REAL *) malloc(k * sizeof(*conv_mem)); + colno_mem = (int *) malloc(k * sizeof(*colno_mem)); + + unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); + RNGType rng(seed); + boost::random::uniform_real_distribution<> urdist1(-1, 1); + + Point p(d); + typename std::vector::iterator pit; + MT V(k, d); + unsigned int j, count_row,it=0; + std::vector indices; + Polytope VP; + VT b = VT::Ones(k); + + for (unsigned int i = 0; i < k; ++i) { + for (int j = 0; j < d; ++j) { + V(i, j) = urdist1(rng); + } + } + if(k==d+1){ + VP.init(d, V, b); + return VP; + } + + MT V2(k,d); + V2 = V; + indices.clear(); + while(it<20) { + V.resize(V2.rows(), d); + V = V2; + for (int i = 0; i < indices.size(); ++i) { + V.conservativeResize(V.rows()+1, d); + for (int j = 0; j < d; ++j) { + V(V.rows()-1, j) = urdist1(rng); + } + } + indices.clear(); + V2.resize(k, d); + V2 = V; + + for (int i = 0; i < k; ++i) { + for (int j = 0; j < d; ++j) { + p.set_coord(j, V(i, j)); + } + removeRow(V2, i); + if (memLP_Vpoly(V2, p, conv_mem, colno_mem)){ + indices.push_back(i); + } + V2.resize(k, d); + V2 = V; + } + if (indices.size()==0) { + VP.init(d, V, b); + return VP; + } + V2.resize(k - indices.size(), d); + count_row =0; + for (int i = 0; i < k; ++i) { + if(std::find(indices.begin(), indices.end(), i) != indices.end()) { + continue; + } else { + for (int j = 0; j < d; ++j) V2(count_row, j) = V(i,j); + count_row++; + } + } + it++; + } + + VP.init(d, V2, VT::Ones(V2.rows())); + free(colno_mem); + free(conv_mem); + + return VP; + +} + +#endif diff --git a/include/generators/z_polytopes_gen.h b/include/generators/z_polytopes_gen.h new file mode 100644 index 000000000..991d8b503 --- /dev/null +++ b/include/generators/z_polytopes_gen.h @@ -0,0 +1,123 @@ +// VolEsti (volume computation and sampling library) + +// Copyright (c) 20012-2018 Vissarion Fisikopoulos +// Copyright (c) 2018 Apostolos Chalkis + +// Licensed under GNU LGPL.3, see LICENCE file + +#ifndef Z_POLYTOPES_GEN_H +#define Z_POLYTOPES_GEN_H + +#include +#include "samplers.h" + +template +Polytope gen_zonotope_gaussian(unsigned int dim, unsigned int m) { + + typedef typename Polytope::MT MT; + typedef typename Polytope::VT VT; + typedef typename Polytope::NT NT; + + unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); + RNGType rng(seed); + boost::normal_distribution<> rdist(0, 1); + boost::normal_distribution<> rdist2(50, 33.3); + + MT A; + VT b; + A.resize(m, dim); + b.resize(m); + Polytope P; + NT rand_gaus; + + for (unsigned int i = 0; i < m; ++i) { + b(i) = 1.0; + for (unsigned int j = 0; j < dim; ++j) { + A(i,j) = rdist(rng); + } + A.row(i)=A.row(i)/A.row(i).norm(); + while(true){ + rand_gaus = rdist2(rng); + if (rand_gaus > 0.0 && rand_gaus<100.0){ + A.row(i) = A.row(i) * rand_gaus; + break; + } + } + } + + P.init(dim, A, b); + return P; +} + + +template +Polytope gen_zonotope_uniform(unsigned int dim, unsigned int m) { + + typedef typename Polytope::MT MT; + typedef typename Polytope::VT VT; + typedef typename Polytope::NT NT; + + unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); + RNGType rng(seed); + boost::normal_distribution<> rdist(0, 1); + boost::random::uniform_real_distribution<> urdist1(0, 100); + + MT A; + VT b; + A.resize(m, dim); + b.resize(m); + Polytope P; + + for (unsigned int i = 0; i < m; ++i) { + b(i) = 1.0; + for (unsigned int j = 0; j < dim; ++j) { + A(i,j) = rdist(rng); + } + A.row(i)=A.row(i)/A.row(i).norm(); + A.row(i) = A.row(i) * urdist1(rng); + } + + P.init(dim, A, b); + return P; +} + + +template +Polytope gen_zonotope_exponential(unsigned int dim, unsigned int m) { + + typedef typename Polytope::MT MT; + typedef typename Polytope::VT VT; + typedef typename Polytope::NT NT; + + unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); + RNGType rng(seed); + boost::normal_distribution<> rdist(0, 1); + boost::normal_distribution<> expdist(1.0/30.0); + + MT A; + VT b; + A.resize(m, dim); + b.resize(m); + Polytope P; + NT rand_exp; + + for (unsigned int i = 0; i < m; ++i) { + b(i) = 1.0; + for (unsigned int j = 0; j < dim; ++j) { + A(i,j) = rdist(rng); + } + A.row(i)=A.row(i)/A.row(i).norm(); + while(true){ + rand_exp = expdist(rng); + if (rand_exp > 0.0 && rand_exp<100.0){ + A.row(i) = A.row(i) * rand_exp; + break; + } + } + } + + P.init(dim, A, b); + return P; +} + +#endif diff --git a/test/VpolyCV_test.cpp b/test/VpolyCV_test.cpp index 87fa8e490..6ee673af3 100644 --- a/test/VpolyCV_test.cpp +++ b/test/VpolyCV_test.cpp @@ -13,7 +13,7 @@ #include "random/normal_distribution.hpp" #include "random/uniform_real_distribution.hpp" #include "volume.h" -#include "polytope_generators.h" +#include "known_polytope_generators.h" #include #include diff --git a/test/VpolyVol_test.cpp b/test/VpolyVol_test.cpp index 8f5ba3c9d..34e8a6c3a 100644 --- a/test/VpolyVol_test.cpp +++ b/test/VpolyVol_test.cpp @@ -13,7 +13,7 @@ #include "random/normal_distribution.hpp" #include "random/uniform_real_distribution.hpp" #include "volume.h" -#include "polytope_generators.h" +#include "known_polytope_generators.h" #include template diff --git a/test/ZonotopeVolCV_test.cpp b/test/ZonotopeVolCV_test.cpp index 0f9eff163..b65604265 100644 --- a/test/ZonotopeVolCV_test.cpp +++ b/test/ZonotopeVolCV_test.cpp @@ -9,7 +9,7 @@ #include #include "Eigen/Eigen" #include "volume.h" -#include "polytope_generators.h" +#include "z_polytopes_gen.h" #include "exact_vols.h" #include @@ -26,7 +26,7 @@ void test_zono_volume(int n, int m, NT tolerance = 0.3) typedef typename Kernel::Point Point; typedef boost::mt19937 RNGType; typedef Zonotope Zonotope; - Zonotope ZP = gen_zonotope(n, m); + Zonotope ZP = gen_zonotope_uniform(n, m); // Setup the parameters int walk_len=1; diff --git a/test/ZonotopeVol_test.cpp b/test/ZonotopeVol_test.cpp index af35577f8..19fce904a 100644 --- a/test/ZonotopeVol_test.cpp +++ b/test/ZonotopeVol_test.cpp @@ -13,7 +13,7 @@ #include "random/normal_distribution.hpp" #include "random/uniform_real_distribution.hpp" #include "volume.h" -#include "polytope_generators.h" +#include "z_polytopes_gen.h" #include "exact_vols.h" #include @@ -30,7 +30,7 @@ void test_zono_volume(int n, int m, NT tolerance = 0.15) typedef typename Kernel::Point Point; typedef boost::mt19937 RNGType; typedef Zonotope Zonotope; - Zonotope ZP = gen_zonotope(n, m); + Zonotope ZP = gen_zonotope_uniform(n, m); // Setup the parameters int walk_len=10 + n/10; diff --git a/test/chebychev_test.cpp b/test/chebychev_test.cpp index 942581f0b..27a03e52d 100644 --- a/test/chebychev_test.cpp +++ b/test/chebychev_test.cpp @@ -15,7 +15,7 @@ #include "random/uniform_real_distribution.hpp" #include "volume.h" #include "misc.h" -#include "polytope_generators.h" +#include "known_polytope_generators.h" template NT factorial(NT n) diff --git a/test/cooling_bodies_bill_test.cpp b/test/cooling_bodies_bill_test.cpp index 2d4331d50..fef5faffe 100644 --- a/test/cooling_bodies_bill_test.cpp +++ b/test/cooling_bodies_bill_test.cpp @@ -29,7 +29,7 @@ #include "cooling_balls.h" #include "cooling_hpoly.h" #include "misc.h" -#include "polytope_generators.h" +#include "known_polytope_generators.h" #include template diff --git a/test/generator.cpp b/test/generator.cpp index 0a7b63c57..f61bbb35c 100644 --- a/test/generator.cpp +++ b/test/generator.cpp @@ -12,15 +12,19 @@ #include "random/uniform_int.hpp" #include "random/normal_distribution.hpp" #include "random/uniform_real_distribution.hpp" -#include "convex_bodies/hpolytope.h" -#include "convex_bodies/vpolytope.h" -#include "convex_bodies/zpolytope.h" -#include "polytope_generators.h" +#include "vpolyoracles.h" +#include "hpolytope.h" +#include "vpolytope.h" +#include "zpolytope.h" +#include "known_polytope_generators.h" +#include "h_polytopes_gen.h" +#include "v_polytopes_gen.h" +#include "z_polytopes_gen.h" #include #include template -void create_txt(MT A, VT b, int kind, bool Vpoly) { +void create_txt(MT A, VT b, int kind, bool Vpoly, bool Zono) { int d = A.cols(), m = A.rows(); std::string bar = "_"; @@ -29,14 +33,46 @@ void create_txt(MT A, VT b, int kind, bool Vpoly) { std::string name; std::ofstream outputFile; - if(Vpoly) { - if (kind == 0) { + if(Zono) { + if (kind == 1) { + std::string poly = "zonotope"; + name = poly + bar + std::to_string(d) + bar + std::to_string(m) + ext; + outputFile.open(name); + outputFile << name << "\n"; + outputFile << "Zonotpe\n"; + } else if (kind == 2) { + std::string poly = "zonotope"; + name = poly + bar + std::to_string(d) + bar + std::to_string(m) + ext; + outputFile.open(name); + outputFile << name << "\n"; + outputFile << "Zonotpe\n"; + } else if (kind == 3) { + std::string poly = "zonotope"; + name = poly + bar + std::to_string(d) + bar + std::to_string(m) + ext; + outputFile.open(name); + outputFile << name << "\n"; + outputFile << "Zonotpe\n"; + } else if (kind == 4) { std::string poly = "zonotope"; name = poly + bar + std::to_string(d) + bar + std::to_string(m) + ext; outputFile.open(name); - outputFile< 0) { - Zonotope ZP = gen_zonotope(d, m); - create_txt(ZP.get_mat(), ZP.get_vec(), kind, true); + Zonotope ZP; + switch (kind) { + case -1: + //std::cout<<"uniform"<(d, m); + break; + case 1: + //std::cout<<"uniform"<(d, m); + break; + case 2: + //std::cout<<"gaussian"<(d, m); + break; + case 3: + //std::cout<<"exponential"<(d, m); + break; + } + //Zonotope ZP = gen_zonotope(d, m); + create_txt(ZP.get_mat(), ZP.get_vec(), kind, false, true); } else { std::cout << "Wrong inputs, try -help" << std::endl; exit(-1); } } else if (Hpoly) { Hpolytope HP; - if (cube) { + if (rand) { + if(m>d+1) { + HP = random_hpoly(d, m); + }else { + std::cout << "the number of facets has to be >=d" << std::endl; + exit(-1); + } + }else if (cube) { HP = gen_cube(d, false); } else if (cross) { HP = gen_cross(d, false); @@ -236,10 +349,21 @@ int main(const int argc, const char** argv) { std::cout << "Wrong inputs, try -help" << std::endl; exit(-1); } - create_txt(HP.get_mat(), HP.get_vec(), kind, false); + create_txt(HP.get_mat(), HP.get_vec(), kind, false, false); } else if (Vpoly) { Vpolytope VP; - if (cube) { + if (rand) { + if (body < 0 || body == 1) { + VP = random_vpoly(d, m); + kind = 5; + } else if (body == 2) { + VP = random_vpoly_incube(d, m); + kind = 4; + } else { + std::cout<<"Wrong inputs, try -help"<(d, true); } else if (cross) { VP = gen_cross(d, true); @@ -255,7 +379,7 @@ int main(const int argc, const char** argv) { std::cout<<"Wrong inputs, try -help"< #include diff --git a/test/volume_test.cpp b/test/volume_test.cpp index d6d21fa23..cae626a3f 100644 --- a/test/volume_test.cpp +++ b/test/volume_test.cpp @@ -15,7 +15,7 @@ #include "random/uniform_real_distribution.hpp" #include "volume.h" #include "misc.h" -#include "polytope_generators.h" +#include "known_polytope_generators.h" #include template From 4ede27b1a1e7b0c4dd12371cb30d26debea85ac0 Mon Sep 17 00:00:00 2001 From: tolischal Date: Mon, 17 Feb 2020 16:36:20 +0200 Subject: [PATCH 05/14] update R random generators --- R-proj/R/RcppExports.R | 4 ++-- R-proj/R/gen_rand_vpoly.R | 16 ++++++++++++---- R-proj/R/gen_rand_zonotope.R | 20 +++++++++++++++----- R-proj/src/poly_gen.cpp | 24 +++++++++++++++++++----- R-proj/src/volume.cpp | 4 ++-- include/generators/z_polytopes_gen.h | 6 +++--- 6 files changed, 53 insertions(+), 21 deletions(-) diff --git a/R-proj/R/RcppExports.R b/R-proj/R/RcppExports.R index 58db215af..7ebdb5e6f 100644 --- a/R-proj/R/RcppExports.R +++ b/R-proj/R/RcppExports.R @@ -122,8 +122,8 @@ inner_ball <- function(P) { #' Do not use this function. #' #' @return A numerical matrix describing the requested polytope -poly_gen <- function(kind_gen, Vpoly_gen, dim_gen, m_gen) { - .Call(`_volesti_poly_gen`, kind_gen, Vpoly_gen, dim_gen, m_gen) +poly_gen <- function(kind_gen, Vpoly_gen, Zono_gen, dim_gen, m_gen) { + .Call(`_volesti_poly_gen`, kind_gen, Vpoly_gen, Zono_gen, dim_gen, m_gen) } #' An internal Rccp function for the random rotation of a convex polytope diff --git a/R-proj/R/gen_rand_vpoly.R b/R-proj/R/gen_rand_vpoly.R index 8b8d7e5da..47001dbb5 100644 --- a/R-proj/R/gen_rand_vpoly.R +++ b/R-proj/R/gen_rand_vpoly.R @@ -3,19 +3,27 @@ #' This function can be used to generate a \eqn{d}-dimensional polytope in V-representation with \eqn{m} vertices. We pick \eqn{m} random points from the boundary of the \eqn{d}-dimensional unit hypersphere as vertices. #' #' @param dimension The dimension of the convex polytope. -#' @param m The number of the vertices. +#' @param n_vertices The number of the vertices. +#' @param generator The body that the generator samples uniformly the vertices from: (a) 'cube' or (b) 'sphere'. #' #' @return A polytope class representing a V-polytope. #' @examples #' # generate a 10-dimensional polytope defined as the convex hull of 25 random vertices #' P = gen_rand_vpoly(10, 25) #' @export -gen_rand_vpoly <- function(dimension, m) { +gen_rand_vpoly <- function(dimension, n_vertices, generator = NULL) { kind_gen = 4 - Vpoly_gen = TRUE - Mat = poly_gen(kind_gen, Vpoly_gen, dimension, m) + if(!is.null(generator)){ + if (generator == 'cube'){ + kind_gen = 5 + } else if (generator != 'sphere') { + stop("Wrong generator!") + } + } + + Mat = poly_gen(kind_gen, TRUE, FALSE, dimension, n_vertices) # first column is the vector b b = Mat[,1] diff --git a/R-proj/R/gen_rand_zonotope.R b/R-proj/R/gen_rand_zonotope.R index e13a079be..f35798be9 100644 --- a/R-proj/R/gen_rand_zonotope.R +++ b/R-proj/R/gen_rand_zonotope.R @@ -3,7 +3,8 @@ #' This function can be used to generate a random \eqn{d}-dimensional zonotope defined by the Minkowski sum of \eqn{m} \eqn{d}-dimensional segments. We consider \eqn{m} random directions in \eqn{R^d} and for each direction we pick a random length in \eqn{[(,\sqrt{d}]} in order to define \eqn{m} segments. #' #' @param dimension The dimension of the zonotope. -#' @param NumGen The number of segments that generate the zonotope. +#' @param n_segments The number of segments that generate the zonotope. +#' @param generator The distribution to pick the length of each segment from \eqn{[0,100]}: (a) 'uniform', (b) 'gaussian' or (c) 'exponential'. #' #' @return A polytope class representing a zonotope. #' @@ -11,12 +12,21 @@ #' # generate a 10-dimensional zonotope defined by the Minkowski sum of 20 segments #' P = gen_rand_zonotope(10, 20) #' @export -gen_rand_zonotope <- function(dimension, NumGen) { +gen_rand_zonotope <- function(dimension, n_segments, generator = NULL) { - kind_gen = 0 - Vpoly_gen = FALSE + kind_gen = 1 - Mat = poly_gen(kind_gen, Vpoly_gen, dimension, NumGen) + if (!is.null(generator)) { + if (generator == 'gaussian') { + kind_gen = 2 + } else if (generator == 'exponential') { + kind_gen = 3 + } else if (generator != 'uniform'){ + stop("Wrong generator!") + } + } + + Mat = poly_gen(kind_gen, FALSE, TRUE, dimension, n_segments) # first column is the vector b b = Mat[,1] diff --git a/R-proj/src/poly_gen.cpp b/R-proj/src/poly_gen.cpp index f9a9e618f..03cea06fb 100644 --- a/R-proj/src/poly_gen.cpp +++ b/R-proj/src/poly_gen.cpp @@ -18,7 +18,10 @@ #include "hpolytope.h" #include "vpolytope.h" #include "zpolytope.h" -#include "polytope_generators.h" +#include "known_polytope_generators.h" +#include "h_polytopes_gen.h" +#include "v_polytopes_gen.h" +#include "z_polytopes_gen.h" #include "extractMatPoly.h" //' An internal Rccp function as a polytope generator @@ -33,7 +36,7 @@ //' //' @return A numerical matrix describing the requested polytope // [[Rcpp::export]] -Rcpp::NumericMatrix poly_gen (int kind_gen, bool Vpoly_gen, int dim_gen, int m_gen) { +Rcpp::NumericMatrix poly_gen (int kind_gen, bool Vpoly_gen, bool Zono_gen, int dim_gen, int m_gen) { typedef double NT; typedef Cartesian Kernel; @@ -43,9 +46,17 @@ Rcpp::NumericMatrix poly_gen (int kind_gen, bool Vpoly_gen, int dim_gen, int m_g typedef VPolytope Vpolytope; typedef Zonotope zonotope; - if (kind_gen == 0) { + if (Zono_gen) { + switch (kind_gen) { - return extractMatPoly(gen_zonotope(dim_gen, m_gen)); + case 1: + return extractMatPoly(gen_zonotope_uniform(dim_gen, m_gen)); + case 2: + return extractMatPoly(gen_zonotope_gaussian(dim_gen, m_gen)); + case 3: + return extractMatPoly(gen_zonotope_exponential(dim_gen, m_gen)); + + } } else if (Vpoly_gen) { switch (kind_gen) { @@ -62,6 +73,9 @@ Rcpp::NumericMatrix poly_gen (int kind_gen, bool Vpoly_gen, int dim_gen, int m_g case 4: return extractMatPoly(random_vpoly(dim_gen, m_gen)); + case 5: + return extractMatPoly(random_vpoly_incube(dim_gen, m_gen)); + } } else { switch (kind_gen) { @@ -87,6 +101,6 @@ Rcpp::NumericMatrix poly_gen (int kind_gen, bool Vpoly_gen, int dim_gen, int m_g } } - return Rcpp::NumericMatrix(0, 0); + throw Rcpp::exception("Wrong inputs!"); } diff --git a/R-proj/src/volume.cpp b/R-proj/src/volume.cpp index b8fda75dc..2041083f2 100644 --- a/R-proj/src/volume.cpp +++ b/R-proj/src/volume.cpp @@ -206,7 +206,7 @@ double volume (Rcpp::Reference P, Rcpp::Nullable walk_length = R_ win_len = 3*n*n+400; } else { billiard = true; - win_len = 150; + win_len = 170; NN = 125; } }else if (Rcpp::as(random_walk).compare(std::string("CDHR")) == 0) { @@ -218,7 +218,7 @@ double volume (Rcpp::Reference P, Rcpp::Nullable walk_length = R_ ball_walk = true; } else if (Rcpp::as(random_walk).compare(std::string("BiW"))==0) { billiard = true; - win_len = 150; + win_len = 170; NN = 125; }else { throw Rcpp::exception("Unknown walk type!"); diff --git a/include/generators/z_polytopes_gen.h b/include/generators/z_polytopes_gen.h index 991d8b503..3aa818df0 100644 --- a/include/generators/z_polytopes_gen.h +++ b/include/generators/z_polytopes_gen.h @@ -12,7 +12,7 @@ #include "samplers.h" template -Polytope gen_zonotope_gaussian(unsigned int dim, unsigned int m) { +Polytope gen_zonotope_gaussian(int dim, int m) { typedef typename Polytope::MT MT; typedef typename Polytope::VT VT; @@ -51,7 +51,7 @@ Polytope gen_zonotope_gaussian(unsigned int dim, unsigned int m) { template -Polytope gen_zonotope_uniform(unsigned int dim, unsigned int m) { +Polytope gen_zonotope_uniform(int dim, int m) { typedef typename Polytope::MT MT; typedef typename Polytope::VT VT; @@ -83,7 +83,7 @@ Polytope gen_zonotope_uniform(unsigned int dim, unsigned int m) { template -Polytope gen_zonotope_exponential(unsigned int dim, unsigned int m) { +Polytope gen_zonotope_exponential(int dim, int m) { typedef typename Polytope::MT MT; typedef typename Polytope::VT VT; From 4d11ee7adde4e16a04ff90a220829345f393a997 Mon Sep 17 00:00:00 2001 From: tolischal Date: Mon, 17 Feb 2020 16:42:02 +0200 Subject: [PATCH 06/14] update Rd files --- R-proj/man/gen_rand_vpoly.Rd | 6 ++++-- R-proj/man/gen_rand_zonotope.Rd | 6 ++++-- R-proj/man/poly_gen.Rd | 2 +- 3 files changed, 9 insertions(+), 5 deletions(-) diff --git a/R-proj/man/gen_rand_vpoly.Rd b/R-proj/man/gen_rand_vpoly.Rd index 14d904e64..a994ce0ef 100644 --- a/R-proj/man/gen_rand_vpoly.Rd +++ b/R-proj/man/gen_rand_vpoly.Rd @@ -4,12 +4,14 @@ \alias{gen_rand_vpoly} \title{Generator function for random V-polytopes} \usage{ -gen_rand_vpoly(dimension, m) +gen_rand_vpoly(dimension, n_vertices, generator = NULL) } \arguments{ \item{dimension}{The dimension of the convex polytope.} -\item{m}{The number of the vertices.} +\item{n_vertices}{The number of the vertices.} + +\item{generator}{The body that the generator samples uniformly the vertices from: (a) 'cube' or (b) 'sphere'.} } \value{ A polytope class representing a V-polytope. diff --git a/R-proj/man/gen_rand_zonotope.Rd b/R-proj/man/gen_rand_zonotope.Rd index 67b3f7f45..8278490d4 100644 --- a/R-proj/man/gen_rand_zonotope.Rd +++ b/R-proj/man/gen_rand_zonotope.Rd @@ -4,12 +4,14 @@ \alias{gen_rand_zonotope} \title{Generator function for zonotopes} \usage{ -gen_rand_zonotope(dimension, NumGen) +gen_rand_zonotope(dimension, n_segments, generator = NULL) } \arguments{ \item{dimension}{The dimension of the zonotope.} -\item{NumGen}{The number of segments that generate the zonotope.} +\item{n_segments}{The number of segments that generate the zonotope.} + +\item{generator}{The distribution to pick the length of each segment from \eqn{[0,100]}: (a) 'uniform', (b) 'gaussian' or (c) 'exponential'.} } \value{ A polytope class representing a zonotope. diff --git a/R-proj/man/poly_gen.Rd b/R-proj/man/poly_gen.Rd index d199997d5..2e7ce5e02 100644 --- a/R-proj/man/poly_gen.Rd +++ b/R-proj/man/poly_gen.Rd @@ -4,7 +4,7 @@ \alias{poly_gen} \title{An internal Rccp function as a polytope generator} \usage{ -poly_gen(kind_gen, Vpoly_gen, dim_gen, m_gen) +poly_gen(kind_gen, Vpoly_gen, Zono_gen, dim_gen, m_gen) } \arguments{ \item{kind_gen}{An integer to declare the type of the polytope.} From efd763f737e01549ee6cbcf492ed9c2bd3416217 Mon Sep 17 00:00:00 2001 From: tolischal Date: Mon, 17 Feb 2020 16:49:24 +0200 Subject: [PATCH 07/14] update R volume interface --- R-proj/src/volume.cpp | 72 +++++++++++++++++++++++++------------------ 1 file changed, 42 insertions(+), 30 deletions(-) diff --git a/R-proj/src/volume.cpp b/R-proj/src/volume.cpp index 2041083f2..e4743a2bc 100644 --- a/R-proj/src/volume.cpp +++ b/R-proj/src/volume.cpp @@ -200,36 +200,6 @@ double volume (Rcpp::Reference P, Rcpp::Nullable walk_length = R_ NT C = 2.0, ratio = 1.0-1.0/(NT(n)), frac = 0.1, e, delta = -1.0, lb = 0.1, ub = 0.15, p = 0.75, rmax = 0.0, alpha = 0.2, diam = -1.0; - if (!random_walk.isNotNull()) { - if ( type == 1 ){ - cdhr = true; - win_len = 3*n*n+400; - } else { - billiard = true; - win_len = 170; - NN = 125; - } - }else if (Rcpp::as(random_walk).compare(std::string("CDHR")) == 0) { - cdhr = true; - win_len = 3*n*n+400; - } else if (Rcpp::as(random_walk).compare(std::string("RDHR")) == 0) { - rdhr = true; - } else if (Rcpp::as(random_walk).compare(std::string("BaW"))==0) { - ball_walk = true; - } else if (Rcpp::as(random_walk).compare(std::string("BiW"))==0) { - billiard = true; - win_len = 170; - NN = 125; - }else { - throw Rcpp::exception("Unknown walk type!"); - } - - if (!rounding.isNotNull() && type == 2){ - round = true; - } else { - round = (!rounding.isNotNull()) ? false : Rcpp::as(rounding); - } - if(!algo.isNotNull()){ if (type == 2 || type == 3) { @@ -263,6 +233,48 @@ double volume (Rcpp::Reference P, Rcpp::Nullable walk_length = R_ throw Rcpp::exception("Unknown method!"); } + if (!random_walk.isNotNull()) { + if ( type == 1 ){ + cdhr = true; + if (CB) win_len = 3*n*n+400; + } else { + if (CB) { + billiard = true; + win_len = 170; + NN = 125; + } else { + rdhr = true; + } + } + }else if (Rcpp::as(random_walk).compare(std::string("CDHR")) == 0) { + cdhr = true; + if (CB) win_len = 3*n*n+400; + } else if (Rcpp::as(random_walk).compare(std::string("RDHR")) == 0) { + rdhr = true; + } else if (Rcpp::as(random_walk).compare(std::string("BaW"))==0) { + ball_walk = true; + } else if (Rcpp::as(random_walk).compare(std::string("BiW"))==0) { + if (CG) { + Rcpp::Rcout << "Billiard walk is not supported for CG algorithm. RDHR is used."<(rounding); + } + if(parameters.isNotNull()) { if (Rcpp::as(parameters).containsElementNamed("BW_rad")) { From 0649d2e4760d9414fe3cf5b127a5dd8adcce2133 Mon Sep 17 00:00:00 2001 From: tolischal Date: Mon, 17 Feb 2020 18:00:18 +0200 Subject: [PATCH 08/14] fix generators in both c++ and R interfaces, improve exact_volume check in c++ interface --- R-proj/R/gen_cross.R | 2 +- R-proj/R/gen_cube.R | 2 +- R-proj/R/gen_prod_simplex.R | 2 +- R-proj/R/gen_rand_hpoly.R | 2 +- R-proj/R/gen_simplex.R | 2 +- R-proj/R/gen_skinny_cube.R | 2 +- R-proj/src/volume.cpp | 20 +++++++++++---- include/volume/cooling_balls.h | 7 ++++-- test/generator.cpp | 45 ++++++++++++---------------------- test/vol.cpp | 6 ++++- 10 files changed, 47 insertions(+), 43 deletions(-) diff --git a/R-proj/R/gen_cross.R b/R-proj/R/gen_cross.R index 026735ae0..6101f5ee1 100644 --- a/R-proj/R/gen_cross.R +++ b/R-proj/R/gen_cross.R @@ -25,7 +25,7 @@ gen_cross <- function(dimension, repr) { stop('Not a known representation.') } - Mat = poly_gen(kind_gen, Vpoly_gen, dimension, m_gen) + Mat = poly_gen(kind_gen, Vpoly_gen, FALSE, dimension, m_gen) # first column is the vector b b = Mat[,1] diff --git a/R-proj/R/gen_cube.R b/R-proj/R/gen_cube.R index cce62489c..5ecf1bcc0 100644 --- a/R-proj/R/gen_cube.R +++ b/R-proj/R/gen_cube.R @@ -25,7 +25,7 @@ gen_cube <- function(dimension, repr) { stop('Not a known representation.') } - Mat = poly_gen(kind_gen, Vpoly_gen, dimension, m_gen) + Mat = poly_gen(kind_gen, Vpoly_gen, FALSE, dimension, m_gen) # first column is the vector b b = Mat[,1] diff --git a/R-proj/R/gen_prod_simplex.R b/R-proj/R/gen_prod_simplex.R index 95e506948..2eda50bd9 100644 --- a/R-proj/R/gen_prod_simplex.R +++ b/R-proj/R/gen_prod_simplex.R @@ -16,7 +16,7 @@ gen_prod_simplex <- function(dimension) { m_gen = 0 Vpoly_gen = FALSE - Mat = poly_gen(kind_gen, Vpoly_gen, dimension, m_gen) + Mat = poly_gen(kind_gen, Vpoly_gen, FALSE, dimension, m_gen) # first column is the vector b b = Mat[,1] diff --git a/R-proj/R/gen_rand_hpoly.R b/R-proj/R/gen_rand_hpoly.R index 0093e09e5..9a5059119 100644 --- a/R-proj/R/gen_rand_hpoly.R +++ b/R-proj/R/gen_rand_hpoly.R @@ -15,7 +15,7 @@ gen_rand_hpoly <- function(dimension, m) { kind_gen = 6 Vpoly_gen = FALSE - Mat = poly_gen(kind_gen, Vpoly_gen, dimension, m) + Mat = poly_gen(kind_gen, Vpoly_gen, FALSE, dimension, m) # first column is the vector b b = Mat[,1] diff --git a/R-proj/R/gen_simplex.R b/R-proj/R/gen_simplex.R index b1f04d917..6d930cdd4 100644 --- a/R-proj/R/gen_simplex.R +++ b/R-proj/R/gen_simplex.R @@ -25,7 +25,7 @@ gen_simplex <- function(dimension, repr) { stop('Not a known representation.') } - Mat = poly_gen(kind_gen, Vpoly_gen, dimension, m_gen) + Mat = poly_gen(kind_gen, Vpoly_gen, FALSE, dimension, m_gen) # first column is the vector b b = Mat[,1] diff --git a/R-proj/R/gen_skinny_cube.R b/R-proj/R/gen_skinny_cube.R index 2c2dbc67f..bcf3dda40 100644 --- a/R-proj/R/gen_skinny_cube.R +++ b/R-proj/R/gen_skinny_cube.R @@ -16,7 +16,7 @@ gen_skinny_cube <- function(dimension) { m_gen = 0 Vpoly_gen = FALSE - Mat = poly_gen(kind_gen, Vpoly_gen, dimension, m_gen) + Mat = poly_gen(kind_gen, Vpoly_gen, FALSE, dimension, m_gen) # first column is the vector b b = Mat[,1] diff --git a/R-proj/src/volume.cpp b/R-proj/src/volume.cpp index e4743a2bc..52bb99cdc 100644 --- a/R-proj/src/volume.cpp +++ b/R-proj/src/volume.cpp @@ -255,11 +255,21 @@ double volume (Rcpp::Reference P, Rcpp::Nullable walk_length = R_ ball_walk = true; } else if (Rcpp::as(random_walk).compare(std::string("BiW"))==0) { if (CG) { - Rcpp::Rcout << "Billiard walk is not supported for CG algorithm. RDHR is used."< BallSet; std::vector ratios; Point c = InnerBall.first; + P.normalize(); if (round) { #ifdef VOLESTI_DEBUG @@ -52,7 +53,10 @@ NT vol_cooling_balls(Polytope &P, UParameters &var, AParameters &var_ban, std::p std::pair res = P.ComputeInnerBall(); c = res.first; radius = res.second; - P.comp_diam(var.diameter, radius); + P.normalize(); + if (var.bill_walk) { + P.comp_diam(var.diameter, radius); + } if (var.ball_walk){ var.delta = 4.0 * radius / NT(n); } @@ -63,7 +67,6 @@ NT vol_cooling_balls(Polytope &P, UParameters &var, AParameters &var_ban, std::p // Move the chebychev center to the origin and apply the same shifting to the polytope VT c_e = Eigen::Map(&c.get_coeffs()[0], c.dimension()); P.shift(c_e); - P.normalize(); if ( !get_sequence_of_polyballs(P, BallSet, ratios, N * nu, nu, lb, ub, radius, alpha, var, rmax) ){ return -1.0; diff --git a/test/generator.cpp b/test/generator.cpp index f61bbb35c..d7b775066 100644 --- a/test/generator.cpp +++ b/test/generator.cpp @@ -169,7 +169,7 @@ int main(const int argc, const char** argv) { bool Hpoly = false, Vpoly = false, Zono = false, cube = false, cross = false, simplex = false, rand = false, prod_simplex = false, skinny_cube = false; - int d = 0, m = 0, kind = -1, body = -1; + int d = 0, m = 0, kind = -1; for(int i=1;i(d, m); break; case 1: - //std::cout<<"uniform"<(d, m); break; case 2: - //std::cout<<"gaussian"<(d, m); break; case 3: - //std::cout<<"exponential"<(d, m); break; } - //Zonotope ZP = gen_zonotope(d, m); create_txt(ZP.get_mat(), ZP.get_vec(), kind, false, true); } else { std::cout << "Wrong inputs, try -help" << std::endl; @@ -332,7 +316,7 @@ int main(const int argc, const char** argv) { if(m>d+1) { HP = random_hpoly(d, m); }else { - std::cout << "the number of facets has to be >=d" << std::endl; + std::cout << "the number of facets has to be >= d+1" << std::endl; exit(-1); } }else if (cube) { @@ -353,14 +337,17 @@ int main(const int argc, const char** argv) { } else if (Vpoly) { Vpolytope VP; if (rand) { - if (body < 0 || body == 1) { - VP = random_vpoly(d, m); - kind = 5; - } else if (body == 2) { - VP = random_vpoly_incube(d, m); - kind = 4; + if (m > 0) { + if (kind == 5) { + VP = random_vpoly(d, m); + } else if (kind == 4) { + VP = random_vpoly_incube(d, m); + } else { + std::cout << "Wrong inputs, try -help" << std::endl; + exit(-1); + } } else { - std::cout<<"Wrong inputs, try -help"<(ZP); std::cout<<"Zonotope's exact volume = "< Date: Mon, 17 Feb 2020 21:15:07 +0200 Subject: [PATCH 09/14] fix bug in hpoly zonotope volume approximation --- include/annealing/hpoly_annealing.h | 3 ++- include/annealing/ratio_estimation.h | 2 +- include/volume/cooling_hpoly.h | 4 ++-- 3 files changed, 5 insertions(+), 4 deletions(-) diff --git a/include/annealing/hpoly_annealing.h b/include/annealing/hpoly_annealing.h index d4c4bd3fd..0a9f929bb 100644 --- a/include/annealing/hpoly_annealing.h +++ b/include/annealing/hpoly_annealing.h @@ -147,7 +147,8 @@ bool get_sequence_of_zonopolys(Zonotope &Z, const HPolytope &HP, std::vector randPoints; Point q(n); diff --git a/include/annealing/ratio_estimation.h b/include/annealing/ratio_estimation.h index a5ac62e77..efe5577c0 100644 --- a/include/annealing/ratio_estimation.h +++ b/include/annealing/ratio_estimation.h @@ -7,7 +7,7 @@ #ifndef RATIO_ESTIMATION_H #define RATIO_ESTIMATION_H -#define MAX_ITER_ESTI 500000 +#define MAX_ITER_ESTI 10000000 template bool check_max_error(const NT &a, const NT &b, const NT &error) { diff --git a/include/volume/cooling_hpoly.h b/include/volume/cooling_hpoly.h index 20372f832..57777bd5a 100644 --- a/include/volume/cooling_hpoly.h +++ b/include/volume/cooling_hpoly.h @@ -54,7 +54,7 @@ NT vol_cooling_hpoly (Zonotope &ZP, UParameters &var, AParameters &var_ban, GPar VT Zs_max(2*m); UParameters var3 = var; var3.cdhr_walk = true; - var3.ball_walk = var3.rdhr_walk = false; + var3.ball_walk = var3.rdhr_walk = var3.bill_walk = false; if ( !get_first_poly(ZP, HP, Zs_max, lb, ub, ratio, var3) ) { return -1.0; } @@ -88,7 +88,7 @@ NT vol_cooling_hpoly (Zonotope &ZP, UParameters &var, AParameters &var_ban, GPar if (!window2) { UParameters var2 = var; var2.cdhr_walk = true; - var2.ball_walk = var2.rdhr_walk = false; + var2.ball_walk = var2.rdhr_walk = var2.bill_walk = false; var2.walk_steps = 10+2*n; vol *= esti_ratio_interval(HP, ZP, ratio, er0, win_len, N*nu, prob, var2); } else { From 8942026b8ddaafa536aa694b34106f6ffaa5d646 Mon Sep 17 00:00:00 2001 From: Tolis Date: Mon, 17 Feb 2020 21:46:36 +0200 Subject: [PATCH 10/14] improve c++ documentation (help command) --- test/vol.cpp | 39 ++++++++++++++++++++++----------------- 1 file changed, 22 insertions(+), 17 deletions(-) diff --git a/test/vol.cpp b/test/vol.cpp index 81b246afb..728b87eae 100644 --- a/test/vol.cpp +++ b/test/vol.cpp @@ -125,27 +125,32 @@ int main(const int argc, const char** argv) "-f3, --file3 [filename_type_G] [epsilon] [walk_length] [threads] [num_of_experiments], for Zonotopes\n"<< "-fle, --filele : counting linear extensions of a poset\n"<< //"-c, --cube [dimension] [epsilon] [walk length] [threads] [num_of_experiments]\n"<< - "--exact : the exact volume\n"<< - "--cube : input polytope is a cube\n"<< - "-exact_zono : compute the exact volume of a zonotope\n"<< + "-exact_vol : the exact volume. Only fo zonotopes\n"<< "-r, --round : enables rounding of the polytope as a preprocess\n"<< "-e, --error epsilon : the goal error of approximation\n"<< - "-w, --walk_len [walk_len] : the random walk length (default 10)\n"<< + "-w, --walk_len [walk_len] : the random walk length. The default value is 1 for CB and CG and 10+d/10 for SOB\n"<< "-exp [#exps] : number of experiments (default 1)\n"<< - "-t, --threads #threads : the number of threads to be used\n"<< - "-ΝΝ : use Nearest Neighbor search to compute the boundary oracles\n"<< - "-birk_sym : use symmetry to compute more random points (only for Birkhoff polytopes)\n"<< - "\n-cg : use the practical CG algo\n"<< + "-cg : use Cooling Gaussians algorithm for volume computation. This is the default choice for H-polytopes in dimension >200\n"<< + "-cb : use Cooling Bodies algorithm for volume computation. This is the default choice for V-polytopes, Zonotopes and H-polytopes in dimension <=200\n"<< + "-sob : use Sequence of Balls algorithm for volume computation\n"<< "-w, --walk_len [walk_len] : the random walk length (default 1)\n"<< - "-rdhr : use random directions HnR, default is coordinate directions HnR\n" - "-e, --error epsilon : the goal error of approximation\n"<< - "-bw : use ball walk for sampling\n"<< + "-rdhr : use random directions HnR\n"<< + "-cdhr : use coordinate directions HnR. This is the default choice for H-polytopes\n"<< + "-biw : use Billiard walk. This is the default choice for V-polytopes and zonotopes\n"<< + "-L : the maximum length of the billiard trajectory\n"<< + "-baw : use ball walk\n"<< "-bwr : the radius of the ball walk (default r*chebychev_radius/sqrt(max(1.0, a_i)*dimension\n"<< - "-Win : the size of the open window for the ratios convergence\n"<< - "-C : a constant for the upper boud of variance/mean^2 in schedule annealing\n" - "-N : the number of points to sample in each step of schedule annealing. Default value N = 500*C + dimension^2/2\n"<< - "-frac : the fraction of the total error to spend in the first gaussian (default frac=0.1)\n"<< - "-ratio : parameter of schedule annealing, larger ratio means larger steps in schedule annealing (default 1-1/dimension)\n"<< + "-e, --error epsilon : the goal error of approximation\n"<< + "-win_len : the size of the open window for the ratios convergence (for CB and CG algorithms)\n"<< + "-C : a constant for the upper boud of variance/mean^2 in schedule annealing. The default values is 2 (for CG algorithm)\n" + "-N : the number of points to sample in each step of schedule annealing. The default value N = 500*C + dimension^2/2 (for CG algorithm)\n"<< + "-frac : the fraction of the total error to spend in the first gaussian. The default frac=0.1 (for CG algo)\n"<< + "-ratio : parameter of schedule annealing, larger ratio means larger steps in schedule annealing. The default 1-1/dimension (for CG algorithm)\n"<< + "-lb : lower bound for the volume ratios in CB algorithm. The default values is 0.1\n"<< + "-ub : upper bound for the volume ratios in CB algorithm. The default values is 0.15\n"<< + "-alpha : the significance level for the statistical test in CB algorithm\n"<< + "-nu : the degrees of freedom of t-student to use in the t-tests in CB algorithm. The default value is 10\n"<< + "-nuN : the degrees of freedom of t-student to use in the t-tests and the number of samples to perform the statistical tests in CB algorithm. The default values is 10 and 125 (when -biw is used) or 120 + d*d/10 (when other random walks are used)\n"<< std::endl; return 0; } @@ -344,7 +349,7 @@ int main(const int argc, const char** argv) nexp = atof(argv[++i]); correct = true; } - if(!strcmp(argv[i],"-diameter")){ + if(!strcmp(argv[i],"-L")){ diameter = atof(argv[++i]); correct = true; } From 653c8064b4e3f02fd105e826a0e6f513c0368176 Mon Sep 17 00:00:00 2001 From: Tolis Date: Mon, 17 Feb 2020 22:33:19 +0200 Subject: [PATCH 11/14] fix c++ tests --- test/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 45755ef91..7a7cc01c7 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -112,7 +112,7 @@ else () add_test(NAME cool_bodies_cube COMMAND cool_bodies_bill_test -tc=cube) add_test(NAME cool_bodies_cross COMMAND cool_bodies_bill_test -tc=cross) add_test(NAME cool_bodies_birkhoff COMMAND cool_bodies_bill_test -tc=birk) - add_test(NAME cool_bodies_prod_simplex COMMAND cool_bodies_bill_test -tc=prod_simplex) + #add_test(NAME cool_bodies_prod_simplex COMMAND cool_bodies_bill_test -tc=prod_simplex) add_test(NAME cool_bodies_simplex COMMAND cool_bodies_bill_test -tc=simplex) add_test(NAME cool_bodies_skinny_cube COMMAND cool_bodies_bill_test -tc=skinny_cube) From 71754f331c7ed695052de089865d3673b064bb1a Mon Sep 17 00:00:00 2001 From: Tolis Date: Mon, 17 Feb 2020 22:59:33 +0200 Subject: [PATCH 12/14] fix c++ tests --- test/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 7a7cc01c7..50c8e0eca 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -79,7 +79,7 @@ else () add_test(NAME volume_cube COMMAND volume_test -tc=cube) add_test(NAME volume_cross COMMAND volume_test -tc=cross) add_test(NAME volume_birkhoff COMMAND volume_test -tc=birk) - add_test(NAME volume_prod_simplex COMMAND volume_test -tc=prod_simplex) + #add_test(NAME volume_prod_simplex COMMAND volume_test -tc=prod_simplex) add_test(NAME volume_simplex COMMAND volume_test -tc=simplex) add_test(NAME volume_skinny_cube COMMAND volume_test -tc=skinny_cube) From 9c4f05c819c80d4e28a85fa1f19515b7f71db90c Mon Sep 17 00:00:00 2001 From: tolischal Date: Wed, 19 Feb 2020 14:53:37 +0200 Subject: [PATCH 13/14] improve cpp generator interface --- test/generator.cpp | 388 ++++++++++++++++++++++----------------------- 1 file changed, 193 insertions(+), 195 deletions(-) diff --git a/test/generator.cpp b/test/generator.cpp index d7b775066..22761cebd 100644 --- a/test/generator.cpp +++ b/test/generator.cpp @@ -27,110 +27,96 @@ template void create_txt(MT A, VT b, int kind, bool Vpoly, bool Zono) { int d = A.cols(), m = A.rows(); - std::string bar = "_"; - std::string ext = ".ext"; - std::string ine = ".ine"; - std::string name; + std::string bar = "_", ext = ".ext", ine = ".ine", name, poly; std::ofstream outputFile; - if(Zono) { - if (kind == 1) { - std::string poly = "zonotope"; - name = poly + bar + std::to_string(d) + bar + std::to_string(m) + ext; - outputFile.open(name); - outputFile << name << "\n"; - outputFile << "Zonotpe\n"; - } else if (kind == 2) { - std::string poly = "zonotope"; - name = poly + bar + std::to_string(d) + bar + std::to_string(m) + ext; - outputFile.open(name); - outputFile << name << "\n"; - outputFile << "Zonotpe\n"; - } else if (kind == 3) { - std::string poly = "zonotope"; - name = poly + bar + std::to_string(d) + bar + std::to_string(m) + ext; - outputFile.open(name); - outputFile << name << "\n"; - outputFile << "Zonotpe\n"; - } else if (kind == 4) { - std::string poly = "zonotope"; - name = poly + bar + std::to_string(d) + bar + std::to_string(m) + ext; - outputFile.open(name); - outputFile << name << "\n"; - outputFile << "Zonotpe\n"; - } - }else if(Vpoly) { - if (kind == 4) { - std::string poly = "rvc"; - name = poly + bar + std::to_string(d) + bar + std::to_string(m) + ext; - outputFile.open(name); - outputFile << name << "\n"; - outputFile << "V-representation\n"; - } else if (kind == 5) { - std::string poly = "rvs"; - name = poly + bar + std::to_string(d) + bar + std::to_string(m) + ext; - outputFile.open(name); - outputFile << name << "\n"; - outputFile << "V-representation\n"; - }else if (kind == 1) { - std::string poly = "cube"; - name = poly + bar + std::to_string(d) + ext; - outputFile.open(name); - outputFile<<"cube_"< 3 || kind <1) std::cout << "Wrong declaration of zonotope's generator, try -help" << std::endl; if (m > 0) { Zonotope ZP; switch (kind) { - case -1: - kind = 1; - ZP = gen_zonotope_uniform(d, m); - break; case 1: ZP = gen_zonotope_uniform(d, m); break; @@ -304,6 +282,10 @@ int main(const int argc, const char** argv) { case 3: ZP = gen_zonotope_exponential(d, m); break; + default: + kind = 1; + ZP = gen_zonotope_uniform(d, m); + break; } create_txt(ZP.get_mat(), ZP.get_vec(), kind, false, true); } else { @@ -312,65 +294,81 @@ int main(const int argc, const char** argv) { } } else if (Hpoly) { Hpolytope HP; - if (rand) { - if(m>d+1) { - HP = random_hpoly(d, m); - }else { - std::cout << "the number of facets has to be >= d+1" << std::endl; + switch (kind) { + case 1: + HP = gen_cube(d, false); + break; + case 2: + HP = gen_cross(d, false); + break; + case 3: + HP = gen_simplex(d, false); + break; + case 4: + HP = gen_prod_simplex(d); + break; + case 5: + HP = gen_skinny_cube(d); + break; + case 6: + if (m > d + 1) { + HP = random_hpoly(d, m); + } else { + std::cout << "The number of facets has to be >= d+1" << std::endl; + exit(-1); + } + break; + default: + std::cout << "Wrong inputs, try -help" << std::endl; exit(-1); - } - }else if (cube) { - HP = gen_cube(d, false); - } else if (cross) { - HP = gen_cross(d, false); - } else if (simplex) { - HP = gen_simplex(d, false); - } else if (prod_simplex) { - HP = gen_prod_simplex(d); - } else if (skinny_cube) { - HP = gen_skinny_cube(d); - } else { - std::cout << "Wrong inputs, try -help" << std::endl; - exit(-1); } create_txt(HP.get_mat(), HP.get_vec(), kind, false, false); } else if (Vpoly) { Vpolytope VP; - if (rand) { - if (m > 0) { - if (kind == 5) { - VP = random_vpoly(d, m); - } else if (kind == 4) { + switch (kind) { + case 1: + VP = gen_cube(d, true); + break; + case 2: + VP = gen_cross(d, true); + break; + case 3: + VP = gen_simplex(d, true); + break; + case 4: + if (!randv) { + std::cout<<"No prod_simplex in V-representation can be generated, try -help"< d) { VP = random_vpoly_incube(d, m); } else { - std::cout << "Wrong inputs, try -help" << std::endl; + std::cout << "The number of vertices has to be >=d+1" << std::endl; exit(-1); } - } else { + break; + case 5: + if (!randv) { + std::cout<<"No skinny_cube in V-representation can be generated, try -help"< d) { + VP = random_vpoly(d, m); + } else { + std::cout << "The number of vertices has to be >=d+1" << std::endl; + exit(-1); + } + break; + default: std::cout << "Wrong inputs, try -help" << std::endl; exit(-1); - } - } else if (cube) { - VP = gen_cube(d, true); - } else if (cross) { - VP = gen_cross(d, true); - } else if (simplex) { - VP = gen_simplex(d, true); - } else if (prod_simplex) { - std::cout<<"No prod_simplex in V-representation can be generated, try -help"< Date: Wed, 19 Feb 2020 15:44:06 +0200 Subject: [PATCH 14/14] improve cpp generator interface --- test/generator.cpp | 64 ++++++++++++---------------------------------- 1 file changed, 17 insertions(+), 47 deletions(-) diff --git a/test/generator.cpp b/test/generator.cpp index 22761cebd..799e9a4eb 100644 --- a/test/generator.cpp +++ b/test/generator.cpp @@ -27,92 +27,62 @@ template void create_txt(MT A, VT b, int kind, bool Vpoly, bool Zono) { int d = A.cols(), m = A.rows(); - std::string bar = "_", ext = ".ext", ine = ".ine", name, poly; + std::string bar = "_", ext = ".ext", ine = ".ine", name; std::ofstream outputFile; if (Zono) { - poly = "zonotope"; - name = poly + bar + std::to_string(d) + bar + std::to_string(m) + ext; + name = "zonotope" + bar + std::to_string(d) + bar + std::to_string(m) + ext; outputFile.open(name); outputFile << name << "\n"; outputFile << "Zonotope\n"; } else if (Vpoly) { switch (kind) { case 1: - poly = "cube"; - name = poly + bar + std::to_string(d) + ext; - outputFile.open(name); - outputFile << "cube_" << d << ".ext\n"; + name = "cube" + bar + std::to_string(d) + ext; break; case 2: - poly = "cross"; - name = poly + bar + std::to_string(d) + ext; - outputFile.open(name); - outputFile << "cross_" << d << ".ext\n"; + name = "cross" + bar + std::to_string(d) + ext; break; case 3: - poly = "simplex"; - name = poly + bar + std::to_string(d) + ext; - outputFile.open(name); - outputFile << "simplex_" << d << ".ext\n"; + name = "simplex" + bar + std::to_string(d) + ext; break; case 4: - poly = "rvc"; - name = poly + bar + std::to_string(d) + bar + std::to_string(m) + ext; - outputFile.open(name); - outputFile << name << "\n"; + name = "rvc" + bar + std::to_string(d) + bar + std::to_string(m) + ext; break; case 5: - poly = "rvs"; - name = poly + bar + std::to_string(d) + bar + std::to_string(m) + ext; - outputFile.open(name); - outputFile << name << "\n"; + name = "rvs" + bar + std::to_string(d) + bar + std::to_string(m) + ext; break; default: return; } + outputFile.open(name); + outputFile << name << "\n"; outputFile << "V-representation\n"; } else { switch (kind) { case 1: - poly = "cube"; - name = poly + bar + std::to_string(d) + ine; - outputFile.open(name); - outputFile << "cube_" << d << ".ine\n"; + name = "cube" + bar + std::to_string(d) + ine; break; case 2: - poly = "cross"; - name = poly + bar + std::to_string(d) + ine; - outputFile.open(name); - outputFile << "cross_" << d << ".ine\n"; + name = "cross" + bar + std::to_string(d) + ine; break; case 3: - poly = "simplex"; - name = poly + bar + std::to_string(d) + ine; - outputFile.open(name); - outputFile << "simplex_" << d << ".ine\n"; + name = "simplex" + bar + std::to_string(d) + ine; break; case 4: - poly = "prod_simplex"; - name = poly + bar + std::to_string(d / 2) + bar + std::to_string(d / 2) + ine; - outputFile.open(name); - outputFile << "prod_simplex_" << d / 2 << "_" << d / 2 << ".ine\n"; + name = "prod_simplex" + bar + std::to_string(d / 2) + bar + std::to_string(d / 2) + ine; break; case 5: - poly = "skinny_cube"; - name = poly + bar + std::to_string(d) + ine; - outputFile.open(name); - outputFile << "skinny_cube_" << d << ".ine\n"; + name = "skinny_cube" + bar + std::to_string(d) + ine; break; case 6: - poly = "rhs"; - name = poly + bar + std::to_string(d) + bar + std::to_string(m) + ine; - outputFile.open(name); - outputFile << "rhs_" << d << "_" << m << ".ine\n"; + name = "rhs" + bar + std::to_string(d) + bar + std::to_string(m) + ine; break; default: return; } + outputFile.open(name); + outputFile << name << "\n"; outputFile << "H-representation\n"; } outputFile << "begin\n";