From d424f3e31817cf33099a903ca42706c16627cb4d Mon Sep 17 00:00:00 2001 From: lucaperju <104904648+lucaperju@users.noreply.github.com> Date: Tue, 30 Jul 2024 21:23:44 +0300 Subject: [PATCH] Support for sparse Eigen matrices in Hpolytope class (#327) * indentation and start of sparse support * sparse matrix support for hpolytope and gabw * unit tests and minor fixes * fix conflict bug * remove debug print and change variable names * fix innerBall bug * has_ball flag, separate A_type and E_type in GABW, improve unit test * rounding and volume_cg unit tests for sparse hpoly * sampling test for sparse polytope * minor changes * fix bug, update set_interior_point function * check that interior point is valid --- .gitignore | 1 + include/convex_bodies/hpolytope.h | 129 ++++++++++++------ include/convex_bodies/orderpolytope.h | 4 +- include/generators/h_polytopes_generator.h | 5 +- .../generators/known_polytope_generators.h | 12 +- include/generators/order_polytope_generator.h | 2 +- include/lp_oracles/solve_lp.h | 4 +- include/random_walks/compute_diameter.hpp | 6 +- .../gaussian_accelerated_billiard_walk.hpp | 71 ++++++---- .../uniform_accelerated_billiard_walk.hpp | 7 +- test/CMakeLists.txt | 4 + test/rounding_test.cpp | 46 ++++++- test/sampling_test.cpp | 97 ++++++++++++- test/volume_cg_hpolytope.cpp | 30 ++++ 14 files changed, 326 insertions(+), 92 deletions(-) diff --git a/.gitignore b/.gitignore index 6976e6147..75e5d999c 100644 --- a/.gitignore +++ b/.gitignore @@ -9,3 +9,4 @@ test/Testing/Temporary/CTestCostData.txt build/ .vscode .DS_Store +.cache \ No newline at end of file diff --git a/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h index b01b6307c..5f37be325 100644 --- a/include/convex_bodies/hpolytope.h +++ b/include/convex_bodies/hpolytope.h @@ -41,16 +41,19 @@ bool is_inner_point_nan_inf(VT const& p) /// This class describes a polytope in H-representation or an H-polytope /// i.e. a polytope defined by a set of linear inequalities /// \tparam Point Point type -template +template +< + typename Point, + typename MT_type = Eigen::Matrix +> class HPolytope { public: typedef Point PointType; typedef typename Point::FT NT; typedef typename std::vector::iterator viterator; - //using RowMatrixXd = Eigen::Matrix; - //typedef RowMatrixXd MT; - typedef Eigen::Matrix MT; + typedef MT_type MT; typedef Eigen::Matrix VT; + typedef Eigen::Matrix DenseMT; private: unsigned int _d; //dimension @@ -58,6 +61,7 @@ class HPolytope { VT b; // vector b, s.t.: Ax<=b std::pair _inner_ball; bool normalized = false; // true if the polytope is normalized + bool has_ball = false; public: //TODO: the default implementation of the Big3 should be ok. Recheck. @@ -68,9 +72,15 @@ class HPolytope { { } + template + HPolytope(unsigned d_, DenseMT const& A_, VT const& b_, typename std::enable_if::value, T>::type* = 0) : + _d{d_}, A{A_.sparseView()}, b{b_} + { + } + // Copy constructor - HPolytope(HPolytope const& p) : - _d{p._d}, A{p.A}, b{p.b}, _inner_ball{p._inner_ball}, normalized{p.normalized} + HPolytope(HPolytope const& p) : + _d{p._d}, A{p.A}, b{p.b}, _inner_ball{p._inner_ball}, normalized{p.normalized}, has_ball{p.has_ball} { } @@ -84,10 +94,10 @@ class HPolytope { for (unsigned int i = 1; i < Pin.size(); i++) { b(i - 1) = Pin[i][0]; for (unsigned int j = 1; j < _d + 1; j++) { - A(i - 1, j - 1) = -Pin[i][j]; + A.coeffRef(i - 1, j - 1) = -Pin[i][j]; } } - _inner_ball.second = -1; + has_ball = false; //_inner_ball = ComputeChebychevBall(A, b); } @@ -100,41 +110,60 @@ class HPolytope { void set_InnerBall(std::pair const& innerball) //const { _inner_ball = innerball; + has_ball = true; } void set_interior_point(Point const& r) { _inner_ball.first = r; + if(!is_in(r)) { + std::cerr << "point outside of polytope was set as interior point" << std::endl; + has_ball = false; + return; + } + if(is_normalized()) { + _inner_ball.second = (b - A * r.getCoefficients()).minCoeff(); + } else { + _inner_ball.second = std::numeric_limits::max(); + for(int i = 0; i < num_of_hyperplanes(); ++i) { + NT dist = (b(i) - A.row(i).dot(r.getCoefficients()) ) / A.row(i).norm(); + if(dist < _inner_ball.second) { + _inner_ball.second = dist; + } + } + } + has_ball = true; } //Compute Chebyshev ball of H-polytope P:= Ax<=b - //Use LpSolve library + //First try using max_inscribed_ball + //Use LpSolve library if it fails std::pair ComputeInnerBall() { normalize(); - #ifndef DISABLE_LPSOLVE - _inner_ball = ComputeChebychevBall(A, b); // use lpsolve library - #else - - if (_inner_ball.second <= NT(0)) { - - NT const tol = 1e-08; - std::tuple inner_ball = max_inscribed_ball(A, b, 5000, tol); - - // check if the solution is feasible - if (is_in(Point(std::get<0>(inner_ball))) == 0 || std::get<1>(inner_ball) < tol/2.0 || - std::isnan(std::get<1>(inner_ball)) || std::isinf(std::get<1>(inner_ball)) || - is_inner_point_nan_inf(std::get<0>(inner_ball))) - { - _inner_ball.second = -1.0; - } else - { - _inner_ball.first = Point(std::get<0>(inner_ball)); - _inner_ball.second = std::get<1>(inner_ball); - } + if (!has_ball) { + + has_ball = true; + NT const tol = 1e-08; + std::tuple inner_ball = max_inscribed_ball(A, b, 5000, tol); + + // check if the solution is feasible + if (is_in(Point(std::get<0>(inner_ball))) == 0 || std::get<1>(inner_ball) < tol/2.0 || + std::isnan(std::get<1>(inner_ball)) || std::isinf(std::get<1>(inner_ball)) || + is_inner_point_nan_inf(std::get<0>(inner_ball))) { + + std::cerr << "Failed to compute max inscribed ball, trying to use lpsolve" << std::endl; + #ifndef DISABLE_LPSOLVE + _inner_ball = ComputeChebychevBall(A, b); // use lpsolve library + #else + std::cerr << "lpsolve is disabled, unable to compute inner ball"; + has_ball = false; + #endif + } else { + _inner_ball.first = Point(std::get<0>(inner_ball)); + _inner_ball.second = std::get<1>(inner_ball); } - #endif - + } return _inner_ball; } @@ -185,6 +214,7 @@ class HPolytope { { A = A2; normalized = false; + has_ball = false; } @@ -192,6 +222,7 @@ class HPolytope { void set_vec(VT const& b2) { b = b2; + has_ball = false; } Point get_mean_of_vertices() const @@ -210,7 +241,7 @@ class HPolytope { std::cout << " " << A.rows() << " " << _d << " double" << std::endl; for (unsigned int i = 0; i < A.rows(); i++) { for (unsigned int j = 0; j < _d; j++) { - std::cout << A(i, j) << " "; + std::cout << A.coeff(i, j) << " "; } std::cout << "<= " << b(i) << std::endl; } @@ -493,7 +524,7 @@ class HPolytope { VT& Ar, VT& Av, NT const& lambda_prev, - MT const& AA, + DenseMT const& AA, update_parameters& params) const { @@ -509,7 +540,7 @@ class HPolytope { if(params.hit_ball) { Av.noalias() += (-2.0 * inner_prev) * (Ar / params.ball_inner_norm); } else { - Av.noalias() += (-2.0 * inner_prev) * AA.col(params.facet_prev); + Av.noalias() += ((-2.0 * inner_prev) * AA.col(params.facet_prev)); } sum_nom.noalias() = b - Ar; @@ -630,12 +661,12 @@ class HPolytope { int m = num_of_hyperplanes(); - lamdas.noalias() += A.col(rand_coord_prev) - * (r_prev[rand_coord_prev] - r[rand_coord_prev]); + lamdas.noalias() += (DenseMT)(A.col(rand_coord_prev) + * (r_prev[rand_coord_prev] - r[rand_coord_prev])); NT* data = lamdas.data(); for (int i = 0; i < m; i++) { - NT a = A(i, rand_coord); + NT a = A.coeff(i, rand_coord); if (a == NT(0)) { //std::cout<<"div0"< + void linear_transformIt(T_type const& T) { - A = A * T; + if constexpr (std::is_same::value) { + A = A * T; + } else { + A = (A * T).sparseView(); + } normalized = false; + has_ball = false; } @@ -838,6 +875,7 @@ class HPolytope { void shift(const VT &c) { b -= A*c; + has_ball = false; } @@ -867,12 +905,15 @@ class HPolytope { void normalize() { + if(normalized) + return; NT row_norm; - for (int i = 0; i < num_of_hyperplanes(); ++i) - { + for (int i = 0; i < A.rows(); ++i) { row_norm = A.row(i).norm(); - A.row(i) = A.row(i) / row_norm; - b(i) = b(i) / row_norm; + if (row_norm != 0.0) { + A.row(i) /= row_norm; + b(i) /= row_norm; + } } normalized = true; } @@ -919,7 +960,7 @@ class HPolytope { } template - NT compute_reflection(Point &v, const Point &, MT const &AE, VT const &AEA, NT const &vEv, update_parameters const ¶ms) const { + NT compute_reflection(Point &v, const Point &, DenseMT const &AE, VT const &AEA, NT const &vEv, update_parameters const ¶ms) const { Point a((-2.0 * params.inner_vi_ak) * A.row(params.facet_prev)); VT x = v.getCoefficients(); diff --git a/include/convex_bodies/orderpolytope.h b/include/convex_bodies/orderpolytope.h index 400d19685..73b4e80f2 100644 --- a/include/convex_bodies/orderpolytope.h +++ b/include/convex_bodies/orderpolytope.h @@ -94,8 +94,8 @@ class OrderPolytope { return _A.sparseView(); } - // return the matrix A - MT get_full_mat() const + // return the matrix A + MT get_dense_mat() const { return _A; } diff --git a/include/generators/h_polytopes_generator.h b/include/generators/h_polytopes_generator.h index 0c9293baa..aec3f1e8f 100644 --- a/include/generators/h_polytopes_generator.h +++ b/include/generators/h_polytopes_generator.h @@ -10,6 +10,7 @@ #include #include +#include #include #include @@ -25,9 +26,9 @@ template Polytope random_hpoly(unsigned int dim, unsigned int m, int seed = std::numeric_limits::signaling_NaN()) { - typedef typename Polytope::MT MT; typedef typename Polytope::VT VT; typedef typename Polytope::NT NT; + typedef typename Eigen::Matrix MT; typedef typename Polytope::PointType Point; int rng_seed = std::chrono::system_clock::now().time_since_epoch().count(); @@ -104,7 +105,7 @@ template Polytope skinny_random_hpoly(unsigned int dim, unsigned int m, const bool pre_rounding = false, const NT eig_ratio = NT(1000.0), int seed = std::numeric_limits::signaling_NaN()) { - typedef typename Polytope::MT MT; + typedef typename Eigen::Matrix MT; typedef typename Polytope::VT VT; typedef typename Polytope::PointType Point; diff --git a/include/generators/known_polytope_generators.h b/include/generators/known_polytope_generators.h index 8d014621d..e88963dee 100644 --- a/include/generators/known_polytope_generators.h +++ b/include/generators/known_polytope_generators.h @@ -19,7 +19,7 @@ template Polytope generate_cube(const unsigned int& dim, const bool& Vpoly,typename Polytope::NT scale=1) { - typedef typename Polytope::MT MT; + typedef typename Eigen::Matrix MT; typedef typename Polytope::VT VT; MT A; VT b; @@ -84,7 +84,7 @@ template Polytope generate_cross(const unsigned int &dim, const bool &Vpoly) { unsigned int m; - typedef typename Polytope::MT MT; + typedef typename Eigen::Matrix MT; typedef typename Polytope::VT VT; MT A; @@ -145,7 +145,7 @@ Polytope generate_cross(const unsigned int &dim, const bool &Vpoly) { /// @tparam Polytope Type of returned polytope template Polytope generate_simplex(const unsigned int &dim, const bool &Vpoly){ - typedef typename Polytope::MT MT; + typedef typename Eigen::Matrix MT; typedef typename Polytope::VT VT; MT A; @@ -198,7 +198,7 @@ Polytope generate_prod_simplex(const unsigned int &dim, bool Vpoly = false){ return Perr; } - typedef typename Polytope::MT MT; + typedef typename Eigen::Matrix MT; typedef typename Polytope::VT VT; MT A; @@ -274,7 +274,7 @@ Polytope generate_skinny_cube(const unsigned int &dim, bool Vpoly = false) { return Perr; } - typedef typename Polytope::MT MT; + typedef typename Eigen::Matrix MT; typedef typename Polytope::VT VT; MT A; @@ -322,7 +322,7 @@ Polytope generate_birkhoff(unsigned int const& n) { unsigned int m = n * n; unsigned int d = n * n - 2 * n + 1; - typedef typename Polytope::MT MT; + typedef typename Eigen::Matrix MT; typedef typename Polytope::VT VT; MT A = MT::Zero(m, d); diff --git a/include/generators/order_polytope_generator.h b/include/generators/order_polytope_generator.h index 3cd907bd9..739169e7e 100644 --- a/include/generators/order_polytope_generator.h +++ b/include/generators/order_polytope_generator.h @@ -63,7 +63,7 @@ Polytope get_orderpoly(Poset const &poset) { if constexpr (std::is_same< Polytope, OrderPolytope >::value ) { return OP; } else if constexpr (std::is_same >::value ){ - Polytope HP(OP.dimension(), OP.get_full_mat(), OP.get_vec()); + Polytope HP(OP.dimension(), OP.get_dense_mat(), OP.get_vec()); return HP; } else { throw "Unable to generate an Order Polytope of requested type"; diff --git a/include/lp_oracles/solve_lp.h b/include/lp_oracles/solve_lp.h index f7013e991..a132ca993 100644 --- a/include/lp_oracles/solve_lp.h +++ b/include/lp_oracles/solve_lp.h @@ -80,8 +80,8 @@ std::pair ComputeChebychevBall(const MT &A, const VT &b){ sum=NT(0); for(j=0; j -struct compute_diameter> +template +struct compute_diameter> { template - static NT compute(HPolytope &P) + static NT compute(HPolytope &P) { return NT(2) * std::sqrt(NT(P.dimension())) * P.InnerBall().second; } diff --git a/include/random_walks/gaussian_accelerated_billiard_walk.hpp b/include/random_walks/gaussian_accelerated_billiard_walk.hpp index ea059305a..d4cd7d9df 100644 --- a/include/random_walks/gaussian_accelerated_billiard_walk.hpp +++ b/include/random_walks/gaussian_accelerated_billiard_walk.hpp @@ -60,20 +60,47 @@ struct GaussianAcceleratedBilliardWalk template < - typename Polytope, - typename RandomNumberGenerator + typename Polytope, + typename RandomNumberGenerator, + typename E_type = typename Polytope::DenseMT > struct Walk { typedef typename Polytope::PointType Point; - typedef typename Polytope::MT MT; + typedef typename Polytope::DenseMT DenseMT; typedef typename Polytope::VT VT; typedef typename Point::FT NT; + void computeLcov(const E_type E) + { + if constexpr (std::is_same >::value) { + Eigen::SimplicialLLT lltofE; + lltofE.compute(E); + if (lltofE.info() != Eigen::Success) { + throw std::runtime_error("First Cholesky decomposition failed for sparse matrix!"); + } + Eigen::SparseMatrix I(E.cols(), E.cols()); + I.setIdentity(); + Eigen::SparseMatrix E_inv = lltofE.solve(I); + Eigen::SimplicialLLT lltofEinv; + lltofEinv.compute(E_inv); + if (lltofE.info() != Eigen::Success) { + throw std::runtime_error("Second Cholesky decomposition failed for sparse matrix!"); + } + _L_cov = lltofEinv.matrixL(); + } else { + Eigen::LLT lltOfE(E.llt().solve(E_type::Identity(E.cols(), E.cols()))); // compute the Cholesky decomposition of inv(E) + if (lltOfE.info() != Eigen::Success) { + throw std::runtime_error("Cholesky decomposition failed for dense matrix!"); + } + _L_cov = lltOfE.matrixL(); + } + } + template Walk(GenericPolytope& P, Point const& p, - MT const& E, // covariance matrix representing the Gaussian distribution + E_type const& E, // covariance matrix representing the Gaussian distribution RandomNumberGenerator &rng) { if(!P.is_normalized()) { @@ -81,13 +108,9 @@ struct GaussianAcceleratedBilliardWalk } _update_parameters = update_parameters(); _L = compute_diameter::template compute(P); - Eigen::LLT lltOfE(E.llt().solve(MT::Identity(E.cols(), E.cols()))); // compute the Cholesky decomposition of inv(E) - if (lltOfE.info() != Eigen::Success) { - throw std::runtime_error("Cholesky decomposition failed!"); - } - _L_cov = lltOfE.matrixL(); + computeLcov(E); _E = E; - _AA.noalias() = P.get_mat() * P.get_mat().transpose(); + _AA.noalias() = (DenseMT)(P.get_mat() * P.get_mat().transpose()); _rho = 1000 * P.dimension(); // upper bound for the number of reflections (experimental) initialize(P, p, rng); } @@ -95,7 +118,7 @@ struct GaussianAcceleratedBilliardWalk template Walk(GenericPolytope& P, Point const& p, - MT const& E, // covariance matrix representing the Gaussian distribution + E_type const& E, // covariance matrix representing the Gaussian distribution RandomNumberGenerator &rng, parameters const& params) { @@ -106,13 +129,9 @@ struct GaussianAcceleratedBilliardWalk _L = params.set_L ? params.m_L : compute_diameter ::template compute(P); - Eigen::LLT lltOfE(E.llt().solve(MT::Identity(E.cols(), E.cols()))); // compute the Cholesky decomposition of inv(E) - if (lltOfE.info() != Eigen::Success) { - throw std::runtime_error("Cholesky decomposition failed!"); - } - _L_cov = lltOfE.matrixL(); + computeLcov(E); _E = E; - _AA.noalias() = P.get_mat() * P.get_mat().transpose(); + _AA.noalias() = (DenseMT)(P.get_mat() * P.get_mat().transpose()); _rho = 1000 * P.dimension(); // upper bound for the number of reflections (experimental) initialize(P, p, rng); } @@ -204,8 +223,14 @@ struct GaussianAcceleratedBilliardWalk _lambdas.setZero(P.num_of_hyperplanes()); _Av.setZero(P.num_of_hyperplanes()); _p = p; - _AE.noalias() = P.get_mat() * _E; - _AEA = _AE.cwiseProduct(P.get_mat()).rowwise().sum(); + _AE.noalias() = (DenseMT)(P.get_mat() * _E); + _AEA = _AE.cwiseProduct((DenseMT)P.get_mat()).rowwise().sum(); + /* + _AEA.resize(P.num_of_hyperplanes()); + for(int i = 0; i < P.num_of_hyperplanes(); ++i) + { + _AEA(i) = _AE.row(i).dot(P.get_mat().row(i)); + }*/ _v = GetDirection::apply(n, rng, false); _v = Point(_L_cov.template triangularView() * _v.getCoefficients()); @@ -257,10 +282,10 @@ struct GaussianAcceleratedBilliardWalk Point _p; Point _v; NT _lambda_prev; - MT _AA; - MT _L_cov; // LL' = inv(E) - MT _AE; - MT _E; + DenseMT _AA; + E_type _L_cov; // LL' = inv(E) + DenseMT _AE; + E_type _E; VT _AEA; unsigned int _rho; update_parameters _update_parameters; diff --git a/include/random_walks/uniform_accelerated_billiard_walk.hpp b/include/random_walks/uniform_accelerated_billiard_walk.hpp index fa9c6fb8f..60ab38263 100644 --- a/include/random_walks/uniform_accelerated_billiard_walk.hpp +++ b/include/random_walks/uniform_accelerated_billiard_walk.hpp @@ -12,6 +12,7 @@ #define RANDOM_WALKS_ACCELERATED_IMPROVED_BILLIARD_WALK_HPP #include "sampling/sphere.hpp" +#include // Billiard walk which accelarates each step for uniform distribution @@ -57,7 +58,7 @@ struct AcceleratedBilliardWalk struct Walk { typedef typename Polytope::PointType Point; - typedef typename Polytope::MT MT; + typedef typename Eigen::Matrix MT; typedef typename Point::FT NT; template @@ -69,7 +70,7 @@ struct AcceleratedBilliardWalk _update_parameters = update_parameters(); _L = compute_diameter ::template compute(P); - _AA.noalias() = P.get_mat() * P.get_mat().transpose(); + _AA.noalias() = (MT)(P.get_mat() * P.get_mat().transpose()); _rho = 1000 * P.dimension(); // upper bound for the number of reflections (experimental) initialize(P, p, rng); } @@ -85,7 +86,7 @@ struct AcceleratedBilliardWalk _L = params.set_L ? params.m_L : compute_diameter ::template compute(P); - _AA.noalias() = P.get_mat() * P.get_mat().transpose(); + _AA.noalias() = (MT)(P.get_mat() * P.get_mat().transpose()); _rho = 1000 * P.dimension(); // upper bound for the number of reflections (experimental) initialize(P, p, rng); } diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 61e499a9e..f57fa007d 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -210,6 +210,7 @@ add_test(NAME test_grdhr COMMAND sampling_test -tc=grdhr) add_test(NAME test_gbaw COMMAND sampling_test -tc=gbaw) add_test(NAME test_ghmc COMMAND sampling_test -tc=ghmc) add_test(NAME test_gabw COMMAND sampling_test -tc=gabw) +add_test(NAME test_sparse COMMAND sampling_test -tc=sparse) add_executable (mmcs_test mmcs_test.cpp $) add_test(NAME test_mmcs COMMAND mmcs_test -tc=mmcs) @@ -263,6 +264,7 @@ add_test(NAME volume_cg_hpolytope_birkhoff COMMAND volume_cg_hpolytope -tc=birk) add_test(NAME volume_cg_hpolytope_prod_simplex COMMAND volume_cg_hpolytope -tc=prod_simplex) add_test(NAME volume_cg_hpolytope_simplex COMMAND volume_cg_hpolytope -tc=simplex) add_test(NAME volume_cg_hpolytope_skinny_cube COMMAND volume_cg_hpolytope -tc=skinny_cube) +add_test(NAME volume_cg_hpolytope_sparse_simplex COMMAND volume_cg_hpolytope -tc=sparse_simplex) add_executable (volume_cg_vpolytope volume_cg_vpolytope.cpp $) add_test(NAME volume_cg_vpolytope_cube COMMAND volume_cg_vpolytope -tc=cube) @@ -305,6 +307,8 @@ add_test(NAME round_volumetric_barrier_test COMMAND rounding_test -tc=round_volumetric_barrier_test) add_test(NAME round_vaidya_barrier_test COMMAND rounding_test -tc=round_vaidya_barrier_test) +add_test(NAME round_sparse_test + COMMAND rounding_test -tc=round_sparse) diff --git a/test/rounding_test.cpp b/test/rounding_test.cpp index 5beae0415..b4739b47e 100644 --- a/test/rounding_test.cpp +++ b/test/rounding_test.cpp @@ -57,7 +57,7 @@ void rounding_min_ellipsoid_test(Polytope &HP, { typedef typename Polytope::PointType Point; typedef typename Point::FT NT; - typedef typename Polytope::MT MT; + typedef typename Polytope::DenseMT MT; typedef typename Polytope::VT VT; int d = HP.dimension(); @@ -102,7 +102,7 @@ void rounding_max_ellipsoid_test(Polytope &HP, { typedef typename Polytope::PointType Point; typedef typename Point::FT NT; - typedef typename Polytope::MT MT; + typedef typename Polytope::DenseMT MT; typedef typename Polytope::VT VT; int d = HP.dimension(); @@ -178,7 +178,7 @@ void rounding_log_barrier_test(Polytope &HP, { typedef typename Polytope::PointType Point; typedef typename Point::FT NT; - typedef typename Polytope::MT MT; + typedef typename Polytope::DenseMT MT; typedef typename Polytope::VT VT; int d = HP.dimension(); @@ -208,7 +208,7 @@ void rounding_volumetric_barrier_test(Polytope &HP, { typedef typename Polytope::PointType Point; typedef typename Point::FT NT; - typedef typename Polytope::MT MT; + typedef typename Polytope::DenseMT MT; typedef typename Polytope::VT VT; int d = HP.dimension(); @@ -238,7 +238,7 @@ void rounding_vaidya_barrier_test(Polytope &HP, { typedef typename Polytope::PointType Point; typedef typename Point::FT NT; - typedef typename Polytope::MT MT; + typedef typename Polytope::DenseMT MT; typedef typename Polytope::VT VT; int d = HP.dimension(); @@ -381,6 +381,39 @@ void call_test_svd() { rounding_svd_test(P, 0, 3070.64, 3188.25, 3140.6, 3200.0); } +template +void call_test_sparse() { + typedef Cartesian Kernel; + typedef typename Kernel::Point Point; + typedef HPolytope > Hpolytope; + Hpolytope P; + + //std::cout << "\n--- Testing SVD rounding of sparse H-skinny_cube5" << std::endl; + //P = generate_skinny_cube(5); + //rounding_svd_test(P, 0, 3070.64, 3188.25, 3140.6, 3200.0); + + std::cout << "\n--- Testing log-barrier rounding of sparse H-skinny_cube5" << std::endl; + P = generate_skinny_cube(5); + rounding_log_barrier_test(P, 0, 3070.64, 3188.25, 3262.77, 3200.0); + + std::cout << "\n--- Testing volumetric barrier rounding of sparse H-skinny_cube5" << std::endl; + P = generate_skinny_cube(5); + rounding_volumetric_barrier_test(P, 0, 3070.64, 3188.25, 3262.77, 3200.0); + + std::cout << "\n--- Testing vaidya barrier rounding of sparse H-skinny_cube5" << std::endl; + P = generate_skinny_cube(5); + rounding_vaidya_barrier_test(P, 0, 3070.64, 3188.25, 3262.77, 3200.0); + + std::cout << "\n--- Testing min ellipsoid rounding of sparse H-skinny_cube10" << std::endl; + P = generate_skinny_cube(10); + rounding_min_ellipsoid_test(P, 0, 122550, 108426, 105003.0, 102400.0); + + std::cout << "\n--- Testing max ellipsoid rounding of sparse H-skinny_cube5" << std::endl; + P = generate_skinny_cube(5); + rounding_max_ellipsoid_test(P, 0, 3070.64, 3188.25, 3262.61, 3200.0); +} + + TEST_CASE("round_min_ellipsoid") { call_test_min_ellipsoid(); @@ -410,3 +443,6 @@ TEST_CASE("round_svd") { call_test_svd(); } +TEST_CASE("round_sparse") { + call_test_sparse(); +} diff --git a/test/sampling_test.cpp b/test/sampling_test.cpp index 887d8489d..a03e46684 100644 --- a/test/sampling_test.cpp +++ b/test/sampling_test.cpp @@ -342,12 +342,103 @@ void call_test_gabw(){ MT samples(d, numpoints); unsigned int jj = 0; - for (typename std::list::iterator rpit = randPoints.begin(); rpit!=randPoints.end(); rpit++, jj++) + for (typename std::list::iterator rpit = randPoints.begin(); rpit != randPoints.end(); rpit++, jj++) samples.col(jj) = (*rpit).getCoefficients(); VT score = univariate_psrf(samples); std::cout << "psrf = " << score.maxCoeff() << std::endl; + CHECK(score.maxCoeff() < 1.1); + RNGType Srng(d); + + typedef Eigen::SparseMatrix SparseMT; + typedef HPolytope SparseHpoly; + std::list Points; + + SparseHpoly SP; + SP = generate_skinny_cube(10); + + + std::cout << "--- Testing Gaussian Accelerated Billiard Walk for Sparse Skinny-H-cube10" << std::endl; + + + p = SP.ComputeInnerBall().first; + typedef typename GaussianAcceleratedBilliardWalk::template Walk + < + SparseHpoly, + RNGType, + SparseMT + > sparsewalk; + typedef MultivariateGaussianRandomPointGenerator SparseRandomPointGenerator; + + + ellipsoid = compute_inscribed_ellipsoid + (SP.get_mat(), SP.get_vec(), p.getCoefficients(), 500, std::pow(10, -6.0), std::pow(10, -4.0)); + + const SparseMT SE = get<0>(ellipsoid).sparseView(); + + SparseRandomPointGenerator::apply(SP, p, SE, numpoints, 1, Points, + push_back_policy, Srng); + + jj = 0; + MT Ssamples(d, numpoints); + for (typename std::list::iterator rpit = Points.begin(); rpit != Points.end(); rpit++, jj++) + Ssamples.col(jj) = (*rpit).getCoefficients(); + + score = univariate_psrf(Ssamples); + std::cout << "psrf = " << score.maxCoeff() << std::endl; + + CHECK(score.maxCoeff() < 1.1); + NT delta = (samples - Ssamples).maxCoeff(); + std::cout << "\ndelta = " << delta << std::endl; + CHECK(delta < 0.00001); + +} + +template +void call_test_sparse(){ + typedef Cartesian Kernel; + typedef typename Kernel::Point Point; + typedef HPolytope> Hpolytope; + typedef typename Hpolytope::MT MT; + typedef typename Hpolytope::DenseMT DenseMT; + typedef Eigen::Matrix VT; + Hpolytope P; + unsigned int d = 10; + + std::cout << "--- Testing ABW for sparse H-cube10" << std::endl; + P = generate_cube(d, false); + P.ComputeInnerBall(); + + DenseMT samples = get_samples(P); + + VT score = univariate_psrf(samples); + std::cout << "psrf = " << score.maxCoeff() << std::endl; + + CHECK(score.maxCoeff() < 1.1); + + + std::cout << "--- Testing CDHR for sparse H-cube10" << std::endl; + P = generate_cube(d, false); + P.ComputeInnerBall(); + + samples = get_samples(P); + + score = univariate_psrf(samples); + std::cout << "psrf = " << score.maxCoeff() << std::endl; + + CHECK(score.maxCoeff() < 1.1); + + + std::cout << "--- Testing RDHR for sparse H-cube10" << std::endl; + P = generate_cube(d, false); + P.ComputeInnerBall(); + + samples = get_samples(P); + + score = univariate_psrf(samples); + std::cout << "psrf = " << score.maxCoeff() << std::endl; + CHECK(score.maxCoeff() < 1.1); } @@ -386,3 +477,7 @@ TEST_CASE("ghmc") { TEST_CASE("gabw") { call_test_gabw(); } + +TEST_CASE("sparse") { + call_test_sparse(); +} diff --git a/test/volume_cg_hpolytope.cpp b/test/volume_cg_hpolytope.cpp index 2e9060e45..e924c4ba3 100644 --- a/test/volume_cg_hpolytope.cpp +++ b/test/volume_cg_hpolytope.cpp @@ -251,6 +251,32 @@ void call_test_skinny_cube() { //test_volume(P, 104857600, 104857600.0); } +template +void call_test_sparse_simplex() { + typedef Cartesian Kernel; + typedef typename Kernel::Point Point; + + typedef HPolytope> Hpolytope; + Hpolytope SP; + + std::cout << "--- Testing volume of sparse H-simplex10" << std::endl; + SP = generate_simplex(10, false); + test_volume(SP, + 2.14048 * std::pow(10,-7), + 2.70598 * std::pow(10,-7), + 2.53893e-07, + 1.0 / factorial(10.0)); + + std::cout << "--- Testing volume of sparse H-simplex20" << std::endl; + SP = generate_simplex(20, false); + test_volume(SP, + 2.00646 * std::pow(10,-21), + 4.16845 * std::pow(10,-19), + 5.0348e-19, + 1.0 / factorial(20.0)); +} + + TEST_CASE("cube") { call_test_cube(); @@ -276,3 +302,7 @@ TEST_CASE("simplex") { TEST_CASE("skinny_cube") { call_test_skinny_cube(); } + +TEST_CASE("sparse_simplex") { + call_test_sparse_simplex(); +}