From d3ef6168bd91283c143af5dbd371739e8710bc8b Mon Sep 17 00:00:00 2001 From: Tolis Chalkis Date: Sat, 2 May 2020 09:51:07 +0300 Subject: [PATCH] Improvements and corrections in Rcpp functions (#4) * correct input list algo in Rcpp volume.cpp * Improvements in rounding and rotating R functions (#76) * add seed in rotate_polytope and update Rd files * update R rounding function to return the linear map * fix rounding in V-poly * overloading rotating() function and remove unused comments * implement zonotope volume computation with cb algorithm and hpoly in MMC * fix all bugs in new cooling_hpoly algorithm and improve test function * all tests for zonotopes with c algortihm passed successfully * implmentation of diameter computation of all convex bodies except Vpoly-intersection * add Vpolyintersection diameter computation and volume approximation * new implmentation of rounding * add rounding in zonotope volume approximation with Hpoly in MMC * dispatching in ratio estimation * update ratio estimation functions * implmemrnt dispatching for diameter computation * improve ratio estimation code structure * update declaration of parameters in random walks * update parameter declaretion for both ball and billiard walks * modify Rcpp expose for volume and sampling C++ functions * modify Rcpp rounding function * update rounding, sampling, volume rcpp dunctions: typedefs and includes * modify zonotope approximation Rcpp function. minor changes to Rcpp volume and sampling functions * fix compiler errors * update d files * fix compiler errors in c++ tests * add seed in rounding and direct_sampling Rcpp functions * add a second volume function to use fixed seed in Rcpp interface * fix rcpp interface --- R-proj/src/rounding.cpp | 13 ++++++----- R-proj/src/sample_points.cpp | 1 + R-proj/src/volume.cpp | 12 +++++----- R-proj/src/zonotope_approximation.cpp | 1 + R-proj/tests/testthat/test_Hvol.R | 22 +++++++++++++------ R-proj/tests/testthat/test_Vvol.R | 16 +++++++------- R-proj/tests/testthat/test_Zvol.R | 16 +++++++------- R-proj/tests/testthat/test_rounding.R | 9 ++++---- .../random_walks/uniform_billiard_walk.hpp | 6 +++++ include/volume/rounding.hpp | 2 +- 10 files changed, 59 insertions(+), 39 deletions(-) diff --git a/R-proj/src/rounding.cpp b/R-proj/src/rounding.cpp index 7650556c2..8ffbb898b 100644 --- a/R-proj/src/rounding.cpp +++ b/R-proj/src/rounding.cpp @@ -15,6 +15,7 @@ #include #include #include +#include "random_walks/random_walks.hpp" #include "volume_sequence_of_balls.hpp" #include "volume_cooling_gaussians.hpp" #include "extractMatPoly.h" @@ -91,10 +92,10 @@ Rcpp::List rounding (Rcpp::Reference P, Rcpp::Nullable seed = R_NilValue switch (type) { case 1: { if (cdhr) { - std::pair , NT> res = round_polytope(HP, InnerBall, walkL, + round_res = round_polytope(HP, InnerBall, walkL, rng); } else { - std::pair , NT> res = round_polytope(HP, InnerBall, walkL, + round_res = round_polytope(HP, InnerBall, walkL, rng); } Mat = extractMatPoly(HP); @@ -102,10 +103,10 @@ Rcpp::List rounding (Rcpp::Reference P, Rcpp::Nullable seed = R_NilValue } case 2: { if (cdhr) { - std::pair , NT> res = round_polytope(VP, InnerBall, walkL, + round_res = round_polytope(VP, InnerBall, walkL, rng); } else { - std::pair , NT> res = round_polytope(VP, InnerBall, walkL, + round_res = round_polytope(VP, InnerBall, walkL, rng); } Mat = extractMatPoly(VP); @@ -113,10 +114,10 @@ Rcpp::List rounding (Rcpp::Reference P, Rcpp::Nullable seed = R_NilValue } case 3: { if (cdhr) { - std::pair , NT> res = round_polytope(ZP, InnerBall, walkL, + round_res = round_polytope(ZP, InnerBall, walkL, rng); } else { - std::pair , NT> res = round_polytope(ZP, InnerBall, walkL, + round_res = round_polytope(ZP, InnerBall, walkL, rng); } Mat = extractMatPoly(ZP); diff --git a/R-proj/src/sample_points.cpp b/R-proj/src/sample_points.cpp index 57d8a5e5b..dec35d8d6 100644 --- a/R-proj/src/sample_points.cpp +++ b/R-proj/src/sample_points.cpp @@ -14,6 +14,7 @@ #include #include #include +#include "random_walks/random_walks.hpp" #include "volume_sequence_of_balls.hpp" #include "volume_cooling_gaussians.hpp" #include "sample_only.h" diff --git a/R-proj/src/volume.cpp b/R-proj/src/volume.cpp index b2909d1f9..a47e97ea7 100644 --- a/R-proj/src/volume.cpp +++ b/R-proj/src/volume.cpp @@ -13,6 +13,8 @@ #include #include #include +#include "random_walks/random_walks.hpp" +#include "volume_cooling_gaussians.hpp" #include "volume_sequence_of_balls.hpp" #include "volume_cooling_gaussians.hpp" #include "volume_cooling_balls.hpp" @@ -47,11 +49,11 @@ double generic_volume(Polytope& P, RNGType &rng, unsigned int walk_length, NT e, NT vol; if (CG) { if (cdhr) { - vol = volume_gaussian_annealing(P, rng, e, walk_length); + vol = volume_cooling_gaussians(P, rng, e, walk_length); } else if (rdhr) { - vol = volume_gaussian_annealing(P, rng, e, walk_length); + vol = volume_cooling_gaussians(P, rng, e, walk_length); } else { - vol = volume_gaussian_annealing(P, rng, e, walk_length); + vol = volume_cooling_gaussians(P, rng, e, walk_length); } } else if (CB) { if (cdhr) { @@ -261,12 +263,12 @@ double volume (Rcpp::Reference P, hpoly = Rcpp::as(Rcpp::as(settings)["hpoly"]); if (hpoly && !CB) Rf_warning("flag 'hpoly' can be used to only in MMC of CB algorithm for zonotopes."); - } else if (ZP.num_of_generators() / ZP.dimension() < 5 ) { + } else if (ZP.num_of_generators() / ZP.dimension() < 5) { hpoly = true; } else { hpoly = false; } - if (hpoly) { + if (hpoly && CB) { if (cdhr) { return volume_cooling_hpoly(ZP, rng, e, walkL, win_len); } else if (rdhr) { diff --git a/R-proj/src/zonotope_approximation.cpp b/R-proj/src/zonotope_approximation.cpp index 1b0989ba7..eabf7c174 100644 --- a/R-proj/src/zonotope_approximation.cpp +++ b/R-proj/src/zonotope_approximation.cpp @@ -11,6 +11,7 @@ #include #include #include +#include "random_walks/random_walks.hpp" #include "volume_sequence_of_balls.hpp" #include "volume_cooling_gaussians.hpp" #include "volume_cooling_balls.hpp" diff --git a/R-proj/tests/testthat/test_Hvol.R b/R-proj/tests/testthat/test_Hvol.R index 3d698421d..03f866290 100644 --- a/R-proj/tests/testthat/test_Hvol.R +++ b/R-proj/tests/testthat/test_Hvol.R @@ -3,14 +3,16 @@ context("H-polytopes' volume test") library(volesti) -Hruntest <- function(P, name_string, exactvol, tol, num_of_exps, alg){ +Hruntest <- function(P, name_string, exactvol, tol, num_of_exps, alg, seed){ vol = 0 for (j in 1:num_of_exps) { if (alg == "CB") { - vol = vol + volume(P) + vol = vol + volume(P, seed = seed) + } else if (alg == "SOB") { + vol = vol + volume(P, settings = list("algorithm" = "SOB"), seed = seed) } else { - vol = vol + volume(P, settings = list("algorithm" = "CG")) + vol = vol + volume(P, settings = list("algorithm" = "CG"), seed = seed) } } vol = vol / num_of_exps @@ -30,7 +32,7 @@ for (i in 1:2) { if (i==1) { algo = 'CG' - num_of_exps = 20 + num_of_exps = 10 } else { algo = 'CB' num_of_exps = 10 @@ -39,20 +41,26 @@ for (i in 1:2) { test_that("Volume H-cube10", { P = gen_cube(10, 'H') - res = Hruntest(P, 'H-cube10', 1024, 0.2, num_of_exps, algo) + res = Hruntest(P, 'H-cube10', 1024, 0.2, num_of_exps, algo, 5) expect_equal(res, 1) }) test_that("Volume H-cross5", { P = gen_cross(5, 'H') - res = Hruntest(P, 'H-cross10', 0.2666667, 0.2, num_of_exps, algo) + res = Hruntest(P, 'H-cross10', 0.2666667, 0.2, num_of_exps, algo, 5) expect_equal(res, 1) }) test_that("Volume H-prod_simplex_5_5", { P = gen_prod_simplex(5) - res = Hruntest(P, 'H-prod_simplex_5_5', (1/prod(1:5))^2, 0.2, num_of_exps, algo) + res = Hruntest(P, 'H-prod_simplex_5_5', (1/prod(1:5))^2, 0.2, num_of_exps, algo, 5) + expect_equal(res, 1) + }) + + test_that("Volume H-cube10", { + P = gen_cube(10, 'H') + res = Hruntest(P, 'H-cube10', 1024, 0.2, 5, "SOB", 5) expect_equal(res, 1) }) diff --git a/R-proj/tests/testthat/test_Vvol.R b/R-proj/tests/testthat/test_Vvol.R index 09964fd64..fcf50aa8e 100644 --- a/R-proj/tests/testthat/test_Vvol.R +++ b/R-proj/tests/testthat/test_Vvol.R @@ -2,14 +2,14 @@ context("V-polytopes' volume test") library(volesti) -Vruntest <- function(P, name_string, exactvol, tol, num_of_exps, algorithm){ +Vruntest <- function(P, name_string, exactvol, tol, num_of_exps, algorithm,seed){ vol = 0 for (j in 1:num_of_exps) { if (algorithm == "CB") { - vol = vol + volume(P) + vol = vol + volume(P, rounding = FALSE, seed = seed) } else { - vol = vol + volume(P, settings = list("algorithm" = "CG", "error" = 0.1)) + vol = vol + volume(P, settings = list("algorithm" = "CG", "error" = 0.1), rounding = FALSE, seed = seed) } } vol = vol / num_of_exps @@ -23,21 +23,21 @@ Vruntest <- function(P, name_string, exactvol, tol, num_of_exps, algorithm){ } cran_only = TRUE -num_of_exps = 6 +num_of_exps = 2 for (i in 1:2) { - + seed = 5 if (i==1) { algo = 'CG' - tol = 0.4 + tol = 0.2 } else { algo = 'CB' - tol = 0.3 + tol = 0.2 } test_that("Volume V-simplex3", { P = gen_simplex(3, 'V') - res = Vruntest(P, 'V-simplex3', 1/prod(1:3), tol, num_of_exps, algo) + res = Vruntest(P, 'V-simplex3', 1/prod(1:3), tol, num_of_exps, algo, seed) expect_equal(res, 1) }) diff --git a/R-proj/tests/testthat/test_Zvol.R b/R-proj/tests/testthat/test_Zvol.R index 3ee3dcd3a..ffda0376c 100644 --- a/R-proj/tests/testthat/test_Zvol.R +++ b/R-proj/tests/testthat/test_Zvol.R @@ -2,15 +2,15 @@ context("Zonotopes' volume test") library(volesti) -Zruntest <- function(P, name_string, tol, num_of_exps, algo){ +Zruntest <- function(P, name_string, tol, num_of_exps, algo, seed){ exactvol = exact_vol(P) vol = 0 for (j in 1:num_of_exps) { if (algo == "CB") { - vol = vol + volume(P, rounding=FALSE) + vol = vol + volume(P, settings = list("hpoly" = FALSE), rounding = FALSE, seed = seed) } else { - vol = vol + volume(P, settings = list("algorithm" = "CG", "error" = 0.1), rounding=TRUE) + vol = vol + volume(P, settings = list("algorithm" = "CG", "error" = 0.1), rounding = FALSE, seed = seed) } } vol = vol / num_of_exps @@ -24,20 +24,20 @@ Zruntest <- function(P, name_string, tol, num_of_exps, algo){ } cran_only = TRUE -num_of_exps = 5 +num_of_exps = 2 for (i in 1:2) { if (i==1) { algo = 'CG' - tol = 0.4 + tol = 0.2 } else { algo = 'CB' - tol = 0.3 + tol = 0.2 } test_that("Volume Zonotope_2_4", { - Z = gen_rand_zonotope(2, 4) - res = Zruntest(Z, 'Zonotope_2_4', tol, num_of_exps, algo) + Z = gen_rand_zonotope(2, 4, seed = 127) + res = Zruntest(Z, 'Zonotope_2_4', tol, num_of_exps, algo, 5) expect_equal(res, 1) }) diff --git a/R-proj/tests/testthat/test_rounding.R b/R-proj/tests/testthat/test_rounding.R index 326a7370d..ad7d03459 100644 --- a/R-proj/tests/testthat/test_rounding.R +++ b/R-proj/tests/testthat/test_rounding.R @@ -2,7 +2,7 @@ context("Rounding test") library(volesti) -testRound <- function(P, exactvol, tol, name_string, num_of_exps, algo, rotation){ +testRound <- function(P, exactvol, tol, name_string, num_of_exps, algo, rotation,seed){ if (rotation) { P = rand_rotate(P) @@ -13,9 +13,9 @@ testRound <- function(P, exactvol, tol, name_string, num_of_exps, algo, rotation vol = 0 for (j in 1:num_of_exps) { if (algo == "CB") { - vol = vol + listHpoly$round_value * volume(listHpoly$P) + vol = vol + listHpoly$round_value * volume(listHpoly$P, seed = seed) } else { - vol = vol + listHpoly$round_value * volume(listHpoly$P, settings=list("algorithm"="CG", "error"=0.1)) + vol = vol + listHpoly$round_value * volume(listHpoly$P, settings=list("algorithm"="CG", "error"=0.1), seed = seed) } } vol = vol / num_of_exps @@ -39,8 +39,9 @@ for (i in 1:2) { test_that("Rounding H-skinny_cube10", { + seed=5 P = gen_skinny_cube(10) - res = testRound(P, 102400, 0.3, 'H-skinny_cube10', num_of_exps, 'CB', FALSE) + res = testRound(P, 102400, 0.3, 'H-skinny_cube10', num_of_exps, 'CB', FALSE,seed) expect_equal(res, 1) }) diff --git a/include/random_walks/uniform_billiard_walk.hpp b/include/random_walks/uniform_billiard_walk.hpp index efeb6025d..80cbdabb6 100644 --- a/include/random_walks/uniform_billiard_walk.hpp +++ b/include/random_walks/uniform_billiard_walk.hpp @@ -384,4 +384,10 @@ private : }; + + + + + + #endif // RANDOM_WALKS_UNIFORM_BILLIARD_WALK_HPP diff --git a/include/volume/rounding.hpp b/include/volume/rounding.hpp index b5d7b7f22..a5ba16c3b 100644 --- a/include/volume/rounding.hpp +++ b/include/volume/rounding.hpp @@ -103,7 +103,7 @@ std::pair< std::pair, NT > round_polytope(Polytope &P, std::pair