From 933fe044618a895a6ee5f771ed4775fd9c5d0d98 Mon Sep 17 00:00:00 2001 From: lucaperju Date: Sat, 20 Jul 2024 19:27:56 +0200 Subject: [PATCH 01/12] indentation and start of sparse support --- include/convex_bodies/hpolytope.h | 10 +++++++--- .../gaussian_accelerated_billiard_walk.hpp | 4 ++-- 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h index b01b6307c..85acc987c 100644 --- a/include/convex_bodies/hpolytope.h +++ b/include/convex_bodies/hpolytope.h @@ -41,7 +41,11 @@ 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 AMat = Eigen::Matrix +> class HPolytope { public: typedef Point PointType; @@ -54,7 +58,7 @@ class HPolytope { private: unsigned int _d; //dimension - MT A; //matrix A + AMat A; //matrix A VT b; // vector b, s.t.: Ax<=b std::pair _inner_ball; bool normalized = false; // true if the polytope is normalized @@ -63,7 +67,7 @@ class HPolytope { //TODO: the default implementation of the Big3 should be ok. Recheck. HPolytope() {} - HPolytope(unsigned d_, MT const& A_, VT const& b_) : + HPolytope(unsigned d_, AMat const& A_, VT const& b_) : _d{d_}, A{A_}, b{b_} { } diff --git a/include/random_walks/gaussian_accelerated_billiard_walk.hpp b/include/random_walks/gaussian_accelerated_billiard_walk.hpp index ea059305a..91deec681 100644 --- a/include/random_walks/gaussian_accelerated_billiard_walk.hpp +++ b/include/random_walks/gaussian_accelerated_billiard_walk.hpp @@ -60,8 +60,8 @@ struct GaussianAcceleratedBilliardWalk template < - typename Polytope, - typename RandomNumberGenerator + typename Polytope, + typename RandomNumberGenerator > struct Walk { From 29ca44efa5f6bc1e2b868a894e0a3d2125536056 Mon Sep 17 00:00:00 2001 From: lucaperju Date: Mon, 22 Jul 2024 16:55:12 +0200 Subject: [PATCH 02/12] sparse matrix support for hpolytope and gabw --- include/convex_bodies/hpolytope.h | 55 ++++++++++++---- 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 | 4 +- include/lp_oracles/solve_lp.h | 4 +- .../gaussian_accelerated_billiard_walk.hpp | 64 +++++++++++++++---- .../uniform_accelerated_billiard_walk.hpp | 13 +++- 8 files changed, 118 insertions(+), 43 deletions(-) diff --git a/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h index 85acc987c..a5131d2d3 100644 --- a/include/convex_bodies/hpolytope.h +++ b/include/convex_bodies/hpolytope.h @@ -44,7 +44,7 @@ bool is_inner_point_nan_inf(VT const& p) template < typename Point, - typename AMat = Eigen::Matrix + typename MT_type = Eigen::Matrix > class HPolytope { public: @@ -53,12 +53,13 @@ class HPolytope { 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 - AMat A; //matrix A + MT A; //matrix A VT b; // vector b, s.t.: Ax<=b std::pair _inner_ball; bool normalized = false; // true if the polytope is normalized @@ -67,13 +68,19 @@ class HPolytope { //TODO: the default implementation of the Big3 should be ok. Recheck. HPolytope() {} - HPolytope(unsigned d_, AMat const& A_, VT const& b_) : + HPolytope(unsigned d_, MT const& A_, VT const& b_) : _d{d_}, A{A_}, b{b_} { } + 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) : + HPolytope(HPolytope const& p) : _d{p._d}, A{p.A}, b{p.b}, _inner_ball{p._inner_ball}, normalized{p.normalized} { } @@ -513,7 +520,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() += (DenseMT)((-2.0 * inner_prev) * AA.col(params.facet_prev)); } sum_nom.noalias() = b - Ar; @@ -830,9 +837,14 @@ class HPolytope { // Apply linear transformation, of square matrix T^{-1}, in H-polytope P:= Ax<=b - void linear_transformIt(MT const& T) + template + 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; } @@ -872,11 +884,28 @@ class HPolytope { void normalize() { NT row_norm; - for (int i = 0; i < num_of_hyperplanes(); ++i) - { - row_norm = A.row(i).norm(); - A.row(i) = A.row(i) / row_norm; - b(i) = b(i) / row_norm; + if constexpr (!std::is_same< MT, Eigen::SparseMatrix >::value ) { + for (int i = 0; i < num_of_hyperplanes(); ++i) + { + row_norm = A.row(i).norm(); + A.row(i) = A.row(i) / row_norm; + b(i) = b(i) / row_norm; + } + } else { + for(int i = 0; i < num_of_hyperplanes(); ++i) + { + row_norm = 0.0; + for(typename Eigen::SparseMatrix::InnerIterator it(A, i); it; ++it) { + row_norm += it.value() * it.value(); + } + row_norm = std::sqrt(row_norm); + if(row_norm != 0.0) { + for (typename Eigen::SparseMatrix::InnerIterator it(A, i); it; ++it) { + it.valueRef() /= row_norm; + } + b(i) = b(i) / row_norm; + } + } } normalized = true; } 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..25f58cc67 100644 --- a/include/generators/order_polytope_generator.h +++ b/include/generators/order_polytope_generator.h @@ -62,8 +62,8 @@ Polytope get_orderpoly(Poset const &poset) { OrderPolytope OP(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()); + } else { + 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 >::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(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 for dense matrix!"); + } + _L_cov = lltOfE.matrixL(); + } + } + template Walk(GenericPolytope& P, Point const& p, @@ -81,13 +107,13 @@ 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(); + if constexpr (std::is_same>::value) { + _AA = P.get_mat() * P.get_mat().transpose(); + } else { + _AA.noalias() = P.get_mat() * P.get_mat().transpose(); + } _rho = 1000 * P.dimension(); // upper bound for the number of reflections (experimental) initialize(P, p, rng); } @@ -106,13 +132,13 @@ 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(); + if constexpr (std::is_same>::value) { + _AA = P.get_mat() * P.get_mat().transpose(); + } else { + _AA.noalias() = 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 +230,18 @@ 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(); + if constexpr (std::is_same>::value) { + _AE = P.get_mat() * _E; + } else { + _AE.noalias() = P.get_mat() * _E; + } + //_AEA = _AE.cwiseProduct(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()); diff --git a/include/random_walks/uniform_accelerated_billiard_walk.hpp b/include/random_walks/uniform_accelerated_billiard_walk.hpp index fa9c6fb8f..f4c7c8ea8 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 @@ -69,7 +70,11 @@ struct AcceleratedBilliardWalk _update_parameters = update_parameters(); _L = compute_diameter ::template compute(P); - _AA.noalias() = P.get_mat() * P.get_mat().transpose(); + if constexpr (std::is_same>::value) { + _AA = P.get_mat() * P.get_mat().transpose(); + } else { + _AA.noalias() = 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 +90,11 @@ struct AcceleratedBilliardWalk _L = params.set_L ? params.m_L : compute_diameter ::template compute(P); - _AA.noalias() = P.get_mat() * P.get_mat().transpose(); + if constexpr (std::is_same>::value) { + _AA = P.get_mat() * P.get_mat().transpose(); + } else { + _AA.noalias() = P.get_mat() * P.get_mat().transpose(); + } _rho = 1000 * P.dimension(); // upper bound for the number of reflections (experimental) initialize(P, p, rng); } From 83a38beb2b1c19cb8be082f072230a1ad9a53e75 Mon Sep 17 00:00:00 2001 From: lucaperju Date: Mon, 22 Jul 2024 22:24:19 +0200 Subject: [PATCH 03/12] unit tests and minor fixes --- include/convex_bodies/hpolytope.h | 71 +++++++++-------------- include/random_walks/compute_diameter.hpp | 7 ++- test/sampling_test.cpp | 42 +++++++++++++- 3 files changed, 72 insertions(+), 48 deletions(-) diff --git a/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h index a5131d2d3..a2cbf5951 100644 --- a/include/convex_bodies/hpolytope.h +++ b/include/convex_bodies/hpolytope.h @@ -120,32 +120,31 @@ class HPolytope { //Compute Chebyshev ball of H-polytope P:= Ax<=b //Use LpSolve library - std::pair ComputeInnerBall() + std::pair ComputeInnerBall(bool const &force = false) { 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))) - { + if (force || _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))) { + + 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"; _inner_ball.second = -1.0; - } else - { - _inner_ball.first = Point(std::get<0>(inner_ball)); - _inner_ball.second = std::get<1>(inner_ball); - } + #endif + } else { + _inner_ball.first = Point(std::get<0>(inner_ball)); + _inner_ball.second = std::get<1>(inner_ball); } - #endif - + } return _inner_ball; } @@ -221,7 +220,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; } @@ -884,27 +883,11 @@ class HPolytope { void normalize() { NT row_norm; - if constexpr (!std::is_same< MT, Eigen::SparseMatrix >::value ) { - for (int i = 0; i < num_of_hyperplanes(); ++i) - { - row_norm = A.row(i).norm(); - A.row(i) = A.row(i) / row_norm; - b(i) = b(i) / row_norm; - } - } else { - for(int i = 0; i < num_of_hyperplanes(); ++i) - { - row_norm = 0.0; - for(typename Eigen::SparseMatrix::InnerIterator it(A, i); it; ++it) { - row_norm += it.value() * it.value(); - } - row_norm = std::sqrt(row_norm); - if(row_norm != 0.0) { - for (typename Eigen::SparseMatrix::InnerIterator it(A, i); it; ++it) { - it.valueRef() /= row_norm; - } - b(i) = b(i) / row_norm; - } + for (int i = 0; i < A.rows(); ++i) { + row_norm = A.row(i).norm(); + if (row_norm != 0.0) { + A.row(i) /= row_norm; + b(i) /= row_norm; } } normalized = true; diff --git a/include/random_walks/compute_diameter.hpp b/include/random_walks/compute_diameter.hpp index b20aa839a..3dc98c7a0 100644 --- a/include/random_walks/compute_diameter.hpp +++ b/include/random_walks/compute_diameter.hpp @@ -30,12 +30,13 @@ struct compute_diameter }; -template -struct compute_diameter> +template +struct compute_diameter> { template - static NT compute(HPolytope &P) + static NT compute(HPolytope &P) { + std::cout << P.InnerBall().second << "<----" << std::endl; return NT(2) * std::sqrt(NT(P.dimension())) * P.InnerBall().second; } }; diff --git a/test/sampling_test.cpp b/test/sampling_test.cpp index 887d8489d..c93d358b1 100644 --- a/test/sampling_test.cpp +++ b/test/sampling_test.cpp @@ -342,13 +342,53 @@ 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); + + + 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 + > 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, rng); + + jj = 0; + MT samples2(d, numpoints); + for (typename std::list::iterator rpit = Points.begin(); rpit != Points.end(); rpit++, jj++) + samples2.col(jj) = (*rpit).getCoefficients(); + + VT score2 = univariate_psrf(samples2); + std::cout << "psrf = " << score2.maxCoeff() << std::endl; + + CHECK(score2.maxCoeff() < 1.1); + } TEST_CASE("dikin") { From 2dea296639828da1b27db62182508d45bf8217d8 Mon Sep 17 00:00:00 2001 From: lucaperju Date: Mon, 22 Jul 2024 22:41:55 +0200 Subject: [PATCH 04/12] fix conflict bug --- include/generators/order_polytope_generator.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/generators/order_polytope_generator.h b/include/generators/order_polytope_generator.h index 25f58cc67..739169e7e 100644 --- a/include/generators/order_polytope_generator.h +++ b/include/generators/order_polytope_generator.h @@ -62,7 +62,7 @@ Polytope get_orderpoly(Poset const &poset) { OrderPolytope OP(poset); if constexpr (std::is_same< Polytope, OrderPolytope >::value ) { return OP; - } else { + } else if constexpr (std::is_same >::value ){ Polytope HP(OP.dimension(), OP.get_dense_mat(), OP.get_vec()); return HP; } else { From 2ad4bfaa56f61dc400b19b177d1132579ec1898d Mon Sep 17 00:00:00 2001 From: lucaperju Date: Mon, 22 Jul 2024 22:45:27 +0200 Subject: [PATCH 05/12] remove debug print and change variable names --- include/random_walks/compute_diameter.hpp | 1 - test/sampling_test.cpp | 9 ++++----- 2 files changed, 4 insertions(+), 6 deletions(-) diff --git a/include/random_walks/compute_diameter.hpp b/include/random_walks/compute_diameter.hpp index 3dc98c7a0..71007d9e4 100644 --- a/include/random_walks/compute_diameter.hpp +++ b/include/random_walks/compute_diameter.hpp @@ -36,7 +36,6 @@ struct compute_diameter> template static NT compute(HPolytope &P) { - std::cout << P.InnerBall().second << "<----" << std::endl; return NT(2) * std::sqrt(NT(P.dimension())) * P.InnerBall().second; } }; diff --git a/test/sampling_test.cpp b/test/sampling_test.cpp index c93d358b1..05284259c 100644 --- a/test/sampling_test.cpp +++ b/test/sampling_test.cpp @@ -380,14 +380,13 @@ void call_test_gabw(){ push_back_policy, rng); jj = 0; - MT samples2(d, numpoints); for (typename std::list::iterator rpit = Points.begin(); rpit != Points.end(); rpit++, jj++) - samples2.col(jj) = (*rpit).getCoefficients(); + samples.col(jj) = (*rpit).getCoefficients(); - VT score2 = univariate_psrf(samples2); - std::cout << "psrf = " << score2.maxCoeff() << std::endl; + score = univariate_psrf(samples); + std::cout << "psrf = " << score.maxCoeff() << std::endl; - CHECK(score2.maxCoeff() < 1.1); + CHECK(score.maxCoeff() < 1.1); } From 44984902630204d1d5640c5325920e228c590c31 Mon Sep 17 00:00:00 2001 From: lucaperju Date: Tue, 23 Jul 2024 20:57:26 +0200 Subject: [PATCH 06/12] fix innerBall bug --- include/convex_bodies/hpolytope.h | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h index a2cbf5951..ab5cff80a 100644 --- a/include/convex_bodies/hpolytope.h +++ b/include/convex_bodies/hpolytope.h @@ -195,6 +195,7 @@ class HPolytope { { A = A2; normalized = false; + _inner_ball.second = -1.0; } @@ -202,6 +203,7 @@ class HPolytope { void set_vec(VT const& b2) { b = b2; + _inner_ball.second = -1.0; } Point get_mean_of_vertices() const @@ -845,6 +847,7 @@ class HPolytope { A = (A * T).sparseView(); } normalized = false; + _inner_ball.second = -1.0; } @@ -853,6 +856,7 @@ class HPolytope { void shift(const VT &c) { b -= A*c; + _inner_ball.second = -1.0; } From bf93e1f8e0f61313648c804438f0ab86c2da5847 Mon Sep 17 00:00:00 2001 From: lucaperju Date: Wed, 24 Jul 2024 07:53:38 +0200 Subject: [PATCH 07/12] has_ball flag, separate A_type and E_type in GABW, improve unit test --- include/convex_bodies/hpolytope.h | 27 +++++++----- .../gaussian_accelerated_billiard_walk.hpp | 44 +++++++++---------- test/sampling_test.cpp | 15 ++++--- 3 files changed, 47 insertions(+), 39 deletions(-) diff --git a/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h index ab5cff80a..bfd506a53 100644 --- a/include/convex_bodies/hpolytope.h +++ b/include/convex_bodies/hpolytope.h @@ -63,6 +63,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. @@ -81,7 +82,7 @@ class HPolytope { // Copy constructor HPolytope(HPolytope const& p) : - _d{p._d}, A{p.A}, b{p.b}, _inner_ball{p._inner_ball}, normalized{p.normalized} + _d{p._d}, A{p.A}, b{p.b}, _inner_ball{p._inner_ball}, normalized{p.normalized}, has_ball{p.has_ball} { } @@ -98,7 +99,7 @@ class HPolytope { A(i - 1, j - 1) = -Pin[i][j]; } } - _inner_ball.second = -1; + has_ball = false; //_inner_ball = ComputeChebychevBall(A, b); } @@ -111,6 +112,7 @@ class HPolytope { void set_InnerBall(std::pair const& innerball) //const { _inner_ball = innerball; + has_ball = true; } void set_interior_point(Point const& r) @@ -120,11 +122,12 @@ class HPolytope { //Compute Chebyshev ball of H-polytope P:= Ax<=b //Use LpSolve library - std::pair ComputeInnerBall(bool const &force = false) + std::pair ComputeInnerBall() { normalize(); - if (force || _inner_ball.second <= NT(0)) { - + if (!has_ball) { + + has_ball = true; NT const tol = 1e-08; std::tuple inner_ball = max_inscribed_ball(A, b, 5000, tol); @@ -138,7 +141,7 @@ class HPolytope { _inner_ball = ComputeChebychevBall(A, b); // use lpsolve library #else std::cerr << "lpsolve is disabled, unable to compute inner ball"; - _inner_ball.second = -1.0; + has_ball = false; #endif } else { _inner_ball.first = Point(std::get<0>(inner_ball)); @@ -195,7 +198,7 @@ class HPolytope { { A = A2; normalized = false; - _inner_ball.second = -1.0; + has_ball = false; } @@ -203,7 +206,7 @@ class HPolytope { void set_vec(VT const& b2) { b = b2; - _inner_ball.second = -1.0; + has_ball = false; } Point get_mean_of_vertices() const @@ -847,7 +850,7 @@ class HPolytope { A = (A * T).sparseView(); } normalized = false; - _inner_ball.second = -1.0; + has_ball = false; } @@ -856,7 +859,7 @@ class HPolytope { void shift(const VT &c) { b -= A*c; - _inner_ball.second = -1.0; + has_ball = false; } @@ -886,6 +889,8 @@ class HPolytope { void normalize() { + if(normalized) + return; NT row_norm; for (int i = 0; i < A.rows(); ++i) { row_norm = A.row(i).norm(); @@ -939,7 +944,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/random_walks/gaussian_accelerated_billiard_walk.hpp b/include/random_walks/gaussian_accelerated_billiard_walk.hpp index ff7a85d60..9efbcf92f 100644 --- a/include/random_walks/gaussian_accelerated_billiard_walk.hpp +++ b/include/random_walks/gaussian_accelerated_billiard_walk.hpp @@ -61,19 +61,21 @@ struct GaussianAcceleratedBilliardWalk template < typename Polytope, - typename RandomNumberGenerator + typename RandomNumberGenerator, + typename E_type = typename Polytope::DenseMT > struct Walk { typedef typename Polytope::PointType Point; - typedef typename Polytope::MT MT; + typedef typename Polytope::MT A_type; + typedef typename Polytope::DenseMT DenseMT; typedef typename Polytope::VT VT; typedef typename Point::FT NT; - void computeLcov(const MT E) + void computeLcov(const E_type E) { - if constexpr (std::is_same >::value) { - Eigen::SimplicialLLT lltofE; + 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!"); @@ -81,14 +83,14 @@ struct GaussianAcceleratedBilliardWalk Eigen::SparseMatrix I(E.cols(), E.cols()); I.setIdentity(); Eigen::SparseMatrix E_inv = lltofE.solve(I); - Eigen::SimplicialLLT lltofEinv; + 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(MT::Identity(E.cols(), E.cols()))); // compute the Cholesky decomposition of inv(E) + 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!"); } @@ -99,7 +101,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) { if(!P.is_normalized()) { @@ -109,7 +111,7 @@ struct GaussianAcceleratedBilliardWalk _L = compute_diameter::template compute(P); computeLcov(E); _E = E; - if constexpr (std::is_same>::value) { + if constexpr (std::is_same>::value) { _AA = P.get_mat() * P.get_mat().transpose(); } else { _AA.noalias() = P.get_mat() * P.get_mat().transpose(); @@ -121,7 +123,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) { @@ -134,7 +136,7 @@ struct GaussianAcceleratedBilliardWalk ::template compute(P); computeLcov(E); _E = E; - if constexpr (std::is_same>::value) { + if constexpr (std::is_same>::value) { _AA = P.get_mat() * P.get_mat().transpose(); } else { _AA.noalias() = P.get_mat() * P.get_mat().transpose(); @@ -230,18 +232,14 @@ struct GaussianAcceleratedBilliardWalk _lambdas.setZero(P.num_of_hyperplanes()); _Av.setZero(P.num_of_hyperplanes()); _p = p; - if constexpr (std::is_same>::value) { - _AE = P.get_mat() * _E; - } else { - _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()); @@ -293,10 +291,10 @@ struct GaussianAcceleratedBilliardWalk Point _p; Point _v; NT _lambda_prev; - MT _AA; - MT _L_cov; // LL' = inv(E) - MT _AE; - MT _E; + A_type _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/test/sampling_test.cpp b/test/sampling_test.cpp index 05284259c..9263a2f82 100644 --- a/test/sampling_test.cpp +++ b/test/sampling_test.cpp @@ -349,7 +349,7 @@ void call_test_gabw(){ std::cout << "psrf = " << score.maxCoeff() << std::endl; CHECK(score.maxCoeff() < 1.1); - + RNGType Srng(d); typedef Eigen::SparseMatrix SparseMT; typedef HPolytope SparseHpoly; @@ -366,7 +366,8 @@ void call_test_gabw(){ typedef typename GaussianAcceleratedBilliardWalk::template Walk < SparseHpoly, - RNGType + RNGType, + SparseMT > sparsewalk; typedef MultivariateGaussianRandomPointGenerator SparseRandomPointGenerator; @@ -377,16 +378,20 @@ void call_test_gabw(){ const SparseMT SE = get<0>(ellipsoid).sparseView(); SparseRandomPointGenerator::apply(SP, p, SE, numpoints, 1, Points, - push_back_policy, rng); + push_back_policy, Srng); jj = 0; + MT Ssamples(d, numpoints); for (typename std::list::iterator rpit = Points.begin(); rpit != Points.end(); rpit++, jj++) - samples.col(jj) = (*rpit).getCoefficients(); + Ssamples.col(jj) = (*rpit).getCoefficients(); - score = univariate_psrf(samples); + 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); } From 08c0f66f91680a78461dbe5420459657f5b9855a Mon Sep 17 00:00:00 2001 From: lucaperju Date: Sat, 27 Jul 2024 15:29:26 +0200 Subject: [PATCH 08/12] rounding and volume_cg unit tests for sparse hpoly --- .gitignore | 1 + include/convex_bodies/hpolytope.h | 8 +++--- test/CMakeLists.txt | 3 ++ test/rounding_test.cpp | 46 +++++++++++++++++++++++++++---- test/volume_cg_hpolytope.cpp | 30 ++++++++++++++++++++ 5 files changed, 79 insertions(+), 9 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 bfd506a53..e30aa25bb 100644 --- a/include/convex_bodies/hpolytope.h +++ b/include/convex_bodies/hpolytope.h @@ -96,7 +96,7 @@ 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]; } } has_ball = false; @@ -645,12 +645,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"<) add_test(NAME volume_cg_vpolytope_cube COMMAND volume_cg_vpolytope -tc=cube) @@ -305,6 +306,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/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(); +} From fe8f92cfc2ece506a41e153293774cd029252c07 Mon Sep 17 00:00:00 2001 From: lucaperju Date: Sat, 27 Jul 2024 20:03:50 +0200 Subject: [PATCH 09/12] sampling test for sparse polytope --- test/CMakeLists.txt | 1 + test/sampling_test.cpp | 51 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 52 insertions(+) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 587966638..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) diff --git a/test/sampling_test.cpp b/test/sampling_test.cpp index 9263a2f82..a03e46684 100644 --- a/test/sampling_test.cpp +++ b/test/sampling_test.cpp @@ -395,6 +395,53 @@ void call_test_gabw(){ } +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); +} + TEST_CASE("dikin") { call_test_dikin(); } @@ -430,3 +477,7 @@ TEST_CASE("ghmc") { TEST_CASE("gabw") { call_test_gabw(); } + +TEST_CASE("sparse") { + call_test_sparse(); +} From 0394b98abda110c0200b0b32487e9d0c6b7c5473 Mon Sep 17 00:00:00 2001 From: lucaperju Date: Tue, 30 Jul 2024 09:10:58 +0200 Subject: [PATCH 10/12] minor changes --- include/convex_bodies/hpolytope.h | 10 +++++----- .../gaussian_accelerated_billiard_walk.hpp | 15 +++------------ .../uniform_accelerated_billiard_walk.hpp | 14 +++----------- 3 files changed, 11 insertions(+), 28 deletions(-) diff --git a/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h index e30aa25bb..6cca06e22 100644 --- a/include/convex_bodies/hpolytope.h +++ b/include/convex_bodies/hpolytope.h @@ -51,8 +51,6 @@ class HPolytope { typedef Point PointType; typedef typename Point::FT NT; typedef typename std::vector::iterator viterator; - //using RowMatrixXd = Eigen::Matrix; - //typedef RowMatrixXd MT; typedef MT_type MT; typedef Eigen::Matrix VT; typedef Eigen::Matrix DenseMT; @@ -118,10 +116,12 @@ class HPolytope { void set_interior_point(Point const& r) { _inner_ball.first = r; + has_ball = false; } //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(); @@ -508,7 +508,7 @@ class HPolytope { VT& Ar, VT& Av, NT const& lambda_prev, - MT const& AA, + DenseMT const& AA, update_parameters& params) const { @@ -524,7 +524,7 @@ class HPolytope { if(params.hit_ball) { Av.noalias() += (-2.0 * inner_prev) * (Ar / params.ball_inner_norm); } else { - Av.noalias() += (DenseMT)((-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; diff --git a/include/random_walks/gaussian_accelerated_billiard_walk.hpp b/include/random_walks/gaussian_accelerated_billiard_walk.hpp index 9efbcf92f..d4cd7d9df 100644 --- a/include/random_walks/gaussian_accelerated_billiard_walk.hpp +++ b/include/random_walks/gaussian_accelerated_billiard_walk.hpp @@ -67,7 +67,6 @@ struct GaussianAcceleratedBilliardWalk struct Walk { typedef typename Polytope::PointType Point; - typedef typename Polytope::MT A_type; typedef typename Polytope::DenseMT DenseMT; typedef typename Polytope::VT VT; typedef typename Point::FT NT; @@ -111,11 +110,7 @@ struct GaussianAcceleratedBilliardWalk _L = compute_diameter::template compute(P); computeLcov(E); _E = E; - if constexpr (std::is_same>::value) { - _AA = P.get_mat() * P.get_mat().transpose(); - } else { - _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); } @@ -136,11 +131,7 @@ struct GaussianAcceleratedBilliardWalk ::template compute(P); computeLcov(E); _E = E; - if constexpr (std::is_same>::value) { - _AA = P.get_mat() * P.get_mat().transpose(); - } else { - _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); } @@ -291,7 +282,7 @@ struct GaussianAcceleratedBilliardWalk Point _p; Point _v; NT _lambda_prev; - A_type _AA; + DenseMT _AA; E_type _L_cov; // LL' = inv(E) DenseMT _AE; E_type _E; diff --git a/include/random_walks/uniform_accelerated_billiard_walk.hpp b/include/random_walks/uniform_accelerated_billiard_walk.hpp index f4c7c8ea8..74185021f 100644 --- a/include/random_walks/uniform_accelerated_billiard_walk.hpp +++ b/include/random_walks/uniform_accelerated_billiard_walk.hpp @@ -58,7 +58,7 @@ struct AcceleratedBilliardWalk struct Walk { typedef typename Polytope::PointType Point; - typedef typename Polytope::MT MT; + typedef typename Polytope::DenseMT MT; typedef typename Point::FT NT; template @@ -70,11 +70,7 @@ struct AcceleratedBilliardWalk _update_parameters = update_parameters(); _L = compute_diameter ::template compute(P); - if constexpr (std::is_same>::value) { - _AA = P.get_mat() * P.get_mat().transpose(); - } else { - _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); } @@ -90,11 +86,7 @@ struct AcceleratedBilliardWalk _L = params.set_L ? params.m_L : compute_diameter ::template compute(P); - if constexpr (std::is_same>::value) { - _AA = P.get_mat() * P.get_mat().transpose(); - } else { - _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); } From 6fe2cf50999b326806fb60ddb3c93d95077b3187 Mon Sep 17 00:00:00 2001 From: lucaperju Date: Tue, 30 Jul 2024 16:05:39 +0200 Subject: [PATCH 11/12] fix bug, update set_interior_point function --- include/convex_bodies/hpolytope.h | 13 ++++++++++++- .../uniform_accelerated_billiard_walk.hpp | 2 +- 2 files changed, 13 insertions(+), 2 deletions(-) diff --git a/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h index 6cca06e22..c60a04ced 100644 --- a/include/convex_bodies/hpolytope.h +++ b/include/convex_bodies/hpolytope.h @@ -116,7 +116,18 @@ class HPolytope { void set_interior_point(Point const& r) { _inner_ball.first = r; - has_ball = false; + 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 diff --git a/include/random_walks/uniform_accelerated_billiard_walk.hpp b/include/random_walks/uniform_accelerated_billiard_walk.hpp index 74185021f..60ab38263 100644 --- a/include/random_walks/uniform_accelerated_billiard_walk.hpp +++ b/include/random_walks/uniform_accelerated_billiard_walk.hpp @@ -58,7 +58,7 @@ struct AcceleratedBilliardWalk struct Walk { typedef typename Polytope::PointType Point; - typedef typename Polytope::DenseMT MT; + typedef typename Eigen::Matrix MT; typedef typename Point::FT NT; template From 52c731ea61f50b79850616d7665b3a78edcccfbc Mon Sep 17 00:00:00 2001 From: lucaperju Date: Tue, 30 Jul 2024 18:17:14 +0200 Subject: [PATCH 12/12] check that interior point is valid --- include/convex_bodies/hpolytope.h | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h index c60a04ced..5f37be325 100644 --- a/include/convex_bodies/hpolytope.h +++ b/include/convex_bodies/hpolytope.h @@ -116,6 +116,11 @@ class HPolytope { 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 {