From d9501a67c0f9674212bcce233f7cb2549ea23ece Mon Sep 17 00:00:00 2001 From: vfisikop <160006479+vfisikop@users.noreply.github.com> Date: Thu, 10 Oct 2024 12:13:03 +0300 Subject: [PATCH] Update volesti source code commit hash ef43e0d18718a158a784b2d5f229decfa4fc3b4d (#28) --- .../convex_bodies/ballintersectconvex.h | 6 + .../correlation_spectrahedron.hpp | 6 + src/volesti/include/convex_bodies/hpolytope.h | 214 ++++++++++++++---- .../include/convex_bodies/orderpolytope.h | 10 + .../spectrahedra/spectrahedron.h | 7 + .../convex_bodies/vpolyintersectvpoly.h | 4 + src/volesti/include/convex_bodies/vpolytope.h | 4 + .../convex_bodies/zonoIntersecthpoly.h | 6 + src/volesti/include/convex_bodies/zpolytope.h | 5 + .../generators/h_polytopes_generator.h | 5 +- .../generators/known_polytope_generators.h | 12 +- .../generators/order_polytope_generator.h | 82 ++++++- src/volesti/include/lp_oracles/solve_lp.h | 4 +- src/volesti/include/misc/poset.h | 65 +++--- .../include/ode_solvers/implicit_midpoint.hpp | 2 + .../include/ode_solvers/oracle_functors.hpp | 73 ++++++ src/volesti/include/preprocess/crhmc/opts.h | 1 + .../estimate_L_smooth_parameter.hpp | 6 +- .../inscribed_ellipsoid_rounding.hpp | 13 +- .../include/preprocess/max_inscribed_ball.hpp | 2 +- .../preprocess/max_inscribed_ellipsoid.hpp | 2 +- .../include/random_walks/compute_diameter.hpp | 6 +- .../additional_units/dynamic_step_size.hpp | 2 +- .../include/random_walks/crhmc/crhmc_walk.hpp | 3 + .../gaussian_accelerated_billiard_walk.hpp | 212 ++++++++++------- .../include/random_walks/gaussian_helpers.hpp | 9 +- .../uniform_accelerated_billiard_walk.hpp | 211 +++++++++++++++-- ...orm_accelerated_billiard_walk_parallel.hpp | 23 +- src/volesti/include/sampling/mmcs.hpp | 2 +- .../include/sampling/parallel_mmcs.hpp | 2 +- .../sampling/random_point_generators.hpp | 18 +- .../random_point_generators_multithread.hpp | 20 +- .../sampling/sample_correlation_matrices.hpp | 15 ++ .../include/volume/volume_cooling_balls.hpp | 6 +- .../volume/volume_cooling_gaussians.hpp | 4 +- 35 files changed, 820 insertions(+), 242 deletions(-) diff --git a/src/volesti/include/convex_bodies/ballintersectconvex.h b/src/volesti/include/convex_bodies/ballintersectconvex.h index b61ae532..992aae1a 100644 --- a/src/volesti/include/convex_bodies/ballintersectconvex.h +++ b/src/volesti/include/convex_bodies/ballintersectconvex.h @@ -50,6 +50,12 @@ class BallIntersectPolytope { return P.get_vec(); } + bool is_normalized () { + return true; + } + + void normalize() {} + int is_in(PointType const& p) const { if (B.is_in(p)==-1) diff --git a/src/volesti/include/convex_bodies/correlation_matrices/correlation_spectrahedron.hpp b/src/volesti/include/convex_bodies/correlation_matrices/correlation_spectrahedron.hpp index a615b33c..92c38dac 100755 --- a/src/volesti/include/convex_bodies/correlation_matrices/correlation_spectrahedron.hpp +++ b/src/volesti/include/convex_bodies/correlation_matrices/correlation_spectrahedron.hpp @@ -261,6 +261,12 @@ class CorrelationSpectrahedron : public Spectrahedron{ MT get_mat() const { return MT::Identity(this->d, this->d); } + + bool is_normalized () { + return true; + } + + void normalize() {} }; #endif //VOLESTI_CONVEX_BODIES_CORRELATION_MATRICES_VOLESTI_CORRELATION_SPECTRAHEDRON_HPP \ No newline at end of file diff --git a/src/volesti/include/convex_bodies/hpolytope.h b/src/volesti/include/convex_bodies/hpolytope.h index d0e894fa..9cabcd9a 100644 --- a/src/volesti/include/convex_bodies/hpolytope.h +++ b/src/volesti/include/convex_bodies/hpolytope.h @@ -6,6 +6,7 @@ //Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018-19 programs. //Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. //Contributed and/or modified by Alexandros Manochis, as part of Google Summer of Code 2020 program. +//Contributed and/or modified by Luca Perju, as part of Google Summer of Code 2024 program. // Licensed under GNU LGPL.3, see LICENCE file @@ -40,22 +41,27 @@ 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 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 + bool has_ball = false; public: //TODO: the default implementation of the Big3 should be ok. Recheck. @@ -66,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} + 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} { } @@ -82,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); } @@ -98,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; } @@ -172,11 +203,18 @@ class HPolytope { return b; } + bool is_normalized () + { + return normalized; + } + // change the matrix A void set_mat(MT const& A2) { A = A2; + normalized = false; + has_ball = false; } @@ -184,6 +222,7 @@ class HPolytope { void set_vec(VT const& b2) { b = b2; + has_ball = false; } Point get_mean_of_vertices() const @@ -202,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; } @@ -485,7 +524,7 @@ class HPolytope { VT& Ar, VT& Av, NT const& lambda_prev, - MT const& AA, + DenseMT const& AA, update_parameters& params) const { @@ -501,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; @@ -528,6 +567,56 @@ class HPolytope { } + + template + std::pair line_positive_intersect(Point const& r, + VT& Ar, + VT& Av, + NT const& lambda_prev, + set_type& distances_set, + AA_type const& AA, + update_parameters& params) const + { + NT inner_prev = params.inner_vi_ak; + NT* Av_data = Av.data(); + + // Av += (-2.0 * inner_prev) * AA.col(params.facet_prev) + + for (Eigen::SparseMatrix::InnerIterator it(AA, params.facet_prev); it; ++it) { + + + // val(row) = (b(row) - Ar(row)) / Av(row) + params.moved_dist + // (val(row) - params.moved_dist) = (b(row) - Ar(row)) / Av(row) + // (val(row) - params.moved_dist) * Av(row) = b(row) - Ar(row) + + *(Av_data + it.row()) += (-2.0 * inner_prev) * it.value(); + + // b(row) - Ar(row) = (old_val(row) - params.moved_dist) * old_Av(row) + // new_val(row) = (b(row) - Ar(row) ) / new_Av(row) + params.moved_dist + // new_val(row) = ((old_val(row) - params.moved_dist) * old_Av(row)) / new_Av(row) + params.moved_dist + + // new_val(row) = (old_val(row) - params.moved_dist) * old_Av(row) / new_Av(row) + params.moved_dist; + // new_val(row) = (old_val(row) - params.moved_dist) * (new_Av(row) + 2.0 * inner_prev * it.value()) / new_Av(row) + params.moved_dist; + // new_val(row) = (old_val(row) - params.moved_dist) * (1 + (2.0 * inner_prev * it.value()) / new_Av(row) ) + params.moved_dist; + + // new_val(row) = old_val(row) + (old_val(row) - params.moved_dist) * 2.0 * inner_prev * it.value() / new_Av(row) + + // val(row) += (val(row) - params.moved_dist) * 2.0 * inner_prev * it.value() / *(Av_data + it.row()); + + NT val = distances_set.get_val(it.row()); + val += (val - params.moved_dist) * 2.0 * inner_prev * it.value() / *(Av_data + it.row()); + distances_set.change_val(it.row(), val, params.moved_dist); + } + + std::pair ans = distances_set.get_min(); + ans.first -= params.moved_dist; + params.inner_vi_ak = *(Av_data + ans.second); + params.facet_prev = ans.second; + + return ans; + } + + template std::pair line_positive_intersect(Point const& r, Point const& v, @@ -541,7 +630,6 @@ class HPolytope { NT lamda = 0; VT sum_nom; - unsigned int j; int m = num_of_hyperplanes(), facet; Ar.noalias() += lambda_prev*Av; @@ -623,12 +711,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; } @@ -830,6 +925,7 @@ class HPolytope { void shift(const VT &c) { b -= A*c; + has_ball = false; } @@ -859,13 +955,17 @@ 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; } void compute_reflection(Point& v, Point const&, int const& facet) const @@ -903,10 +1003,34 @@ class HPolytope { } template - void compute_reflection(Point &v, const Point &, update_parameters const& params) const { + void compute_reflection(Point &v, Point const&, update_parameters const& params) const { + Point a((-2.0 * params.inner_vi_ak) * A.row(params.facet_prev)); + v += a; + } + + // updates the velocity vector v and the position vector p after a reflection + // the real value of p is given by p + moved_dist * v + template + auto compute_reflection(Point &v, Point &p, update_parameters const& params) const + -> std::enable_if_t> && !std::is_same_v, void> { // MT must be in RowMajor format + NT* v_data = v.pointerToData(); + NT* p_data = p.pointerToData(); + for(Eigen::SparseMatrix::InnerIterator it(A, params.facet_prev); it; ++it) { + *(v_data + it.col()) += (-2.0 * params.inner_vi_ak) * it.value(); + *(p_data + it.col()) -= (-2.0 * params.inner_vi_ak * params.moved_dist) * it.value(); + } + } + + template + 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(); + NT new_vEv = vEv - (4.0 * params.inner_vi_ak) * (AE.row(params.facet_prev).dot(x) - params.inner_vi_ak * AEA(params.facet_prev)); v += a; + NT coeff = std::sqrt(vEv / new_vEv); + v = v * coeff; + return coeff; } void update_position_internal(NT&){} diff --git a/src/volesti/include/convex_bodies/orderpolytope.h b/src/volesti/include/convex_bodies/orderpolytope.h index 79f41e16..73b4e80f 100644 --- a/src/volesti/include/convex_bodies/orderpolytope.h +++ b/src/volesti/include/convex_bodies/orderpolytope.h @@ -94,12 +94,22 @@ class OrderPolytope { return _A.sparseView(); } + // return the matrix A + MT get_dense_mat() const + { + return _A; + } + VT get_vec() const { return b; } + bool is_normalized () + { + return _normalized; + } // print polytope in Ax <= b format void print() const diff --git a/src/volesti/include/convex_bodies/spectrahedra/spectrahedron.h b/src/volesti/include/convex_bodies/spectrahedra/spectrahedron.h index e8f34138..37b997ea 100644 --- a/src/volesti/include/convex_bodies/spectrahedra/spectrahedron.h +++ b/src/volesti/include/convex_bodies/spectrahedra/spectrahedron.h @@ -341,6 +341,13 @@ class Spectrahedron { { return MT::Identity(lmi.dimension(), lmi.dimension()); } + + bool is_normalized () + { + return true; + } + + void normalize() {} void resetFlags() { diff --git a/src/volesti/include/convex_bodies/vpolyintersectvpoly.h b/src/volesti/include/convex_bodies/vpolyintersectvpoly.h index a6dd97ca..ec2a7c5d 100644 --- a/src/volesti/include/convex_bodies/vpolyintersectvpoly.h +++ b/src/volesti/include/convex_bodies/vpolyintersectvpoly.h @@ -92,6 +92,10 @@ class IntersectionOfVpoly { return P1.get_vec(); } + bool is_normalized () { + return true; + } + MT get_T() const { return P1.get_mat(); } diff --git a/src/volesti/include/convex_bodies/vpolytope.h b/src/volesti/include/convex_bodies/vpolytope.h index 82ca6689..9100756e 100644 --- a/src/volesti/include/convex_bodies/vpolytope.h +++ b/src/volesti/include/convex_bodies/vpolytope.h @@ -213,6 +213,10 @@ class VPolytope { return b; } + bool is_normalized () { + return true; + } + // change the matrix V void set_mat(const MT &V2) { V = V2; diff --git a/src/volesti/include/convex_bodies/zonoIntersecthpoly.h b/src/volesti/include/convex_bodies/zonoIntersecthpoly.h index bbad43a6..e356db8d 100644 --- a/src/volesti/include/convex_bodies/zonoIntersecthpoly.h +++ b/src/volesti/include/convex_bodies/zonoIntersecthpoly.h @@ -65,6 +65,12 @@ class ZonoIntersectHPoly { return HP.get_vec(); } + bool is_normalized () { + return true; + } + + void normalize() {} + std::pair InnerBall() const { return Z.InnerBall(); diff --git a/src/volesti/include/convex_bodies/zpolytope.h b/src/volesti/include/convex_bodies/zpolytope.h index 1ad8868e..68ce0cf0 100644 --- a/src/volesti/include/convex_bodies/zpolytope.h +++ b/src/volesti/include/convex_bodies/zpolytope.h @@ -299,6 +299,11 @@ class Zonotope { return b; } + bool is_normalized () + { + return true; + } + // change the matrix V void set_mat(MT const& V2) { diff --git a/src/volesti/include/generators/h_polytopes_generator.h b/src/volesti/include/generators/h_polytopes_generator.h index 0c9293ba..aec3f1e8 100644 --- a/src/volesti/include/generators/h_polytopes_generator.h +++ b/src/volesti/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/src/volesti/include/generators/known_polytope_generators.h b/src/volesti/include/generators/known_polytope_generators.h index 8d014621..e88963de 100644 --- a/src/volesti/include/generators/known_polytope_generators.h +++ b/src/volesti/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/src/volesti/include/generators/order_polytope_generator.h b/src/volesti/include/generators/order_polytope_generator.h index bfece556..739169e7 100644 --- a/src/volesti/include/generators/order_polytope_generator.h +++ b/src/volesti/include/generators/order_polytope_generator.h @@ -9,11 +9,26 @@ #ifndef ORDER_POLYTOPES_GEN_H #define ORDER_POLYTOPES_GEN_H +#include #include +#include #include -#include "misc.h" +#include +#include + +#include "misc/misc.h" #include "misc/poset.h" +#include +#include +#include +#include + +#include "generators/boost_random_number_generator.hpp" + +#include "convex_bodies/orderpolytope.h" +#include "convex_bodies/hpolytope.h" + // Instances taken from: https://github.com/ttalvitie/le-counting-practice static const std::unordered_map instances = @@ -38,6 +53,23 @@ static const std::unordered_map instances = }; +// generates a Polytope from a poset +/// @tparam Polytope Type of returned polytope +template +Polytope get_orderpoly(Poset const &poset) { + typedef typename Polytope::PointType Point; + + 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_dense_mat(), OP.get_vec()); + return HP; + } else { + throw "Unable to generate an Order Polytope of requested type"; + } +} + // generates an Order Polytope from an instance name // Instances taken from: https://github.com/ttalvitie/le-counting-practice /// @tparam Polytope Type of returned polytope @@ -45,7 +77,7 @@ template Polytope generate_orderpoly(std::string& instance_name) { std::stringstream in_ss(instances.at(instance_name)); Poset poset = read_poset_from_file_adj_matrix(in_ss).second; - return Polytope(poset); + return get_orderpoly(poset); } // Generates a cube as an Order Polytope @@ -56,8 +88,50 @@ Polytope generate_cube_orderpoly(unsigned int dim) { RV order_relations; Poset poset(dim, order_relations); - Polytope OP(poset); - return OP; + return get_orderpoly(poset); +} + +// Generates a random Order Polytope with given dimension and number of facets +/// @tparam Polytope Type of returned polytope +/// @tparam RNGType RNGType Type +template +Polytope random_orderpoly(unsigned int dim, unsigned int m, int seed = std::numeric_limits::signaling_NaN()) { + + typedef typename Poset::RV RV; + + int rng_seed = std::chrono::system_clock::now().time_since_epoch().count(); + if (!isnan(seed)) { + rng_seed = seed; + } + + typedef BoostRandomNumberGenerator RNG; + RNG rng(dim); + rng.set_seed(rng_seed); + + + std::vector order(dim); + for(int i = 0; i < dim; ++i) { + order[i] = i; + } + boost::mt19937 shuffle_rng(rng_seed); + std::shuffle(order.begin(), order.end(), shuffle_rng); + + + RV order_relations; + for(int i = 0; i < m - 2 * dim; ++i) { + unsigned int x = rng.sample_uidist(); + unsigned int y = rng.sample_uidist(); + while(x == y) { + y = rng.sample_uidist(); + } + if(x > y) + std::swap(x, y); + order_relations.push_back(std::make_pair(order[x], order[y])); + } + + + Poset poset(dim, order_relations); + return get_orderpoly(poset); } #endif diff --git a/src/volesti/include/lp_oracles/solve_lp.h b/src/volesti/include/lp_oracles/solve_lp.h index f7013e99..a132ca99 100644 --- a/src/volesti/include/lp_oracles/solve_lp.h +++ b/src/volesti/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 &res) + { + std::vector > adjList(n); + std::vector indegree(n, 0); + + for(auto x: relations) { + adjList[x.first].push_back(x.second); + indegree[x.second] += 1; + } + + std::queue q; + for(unsigned int i=0; i order; + sorted_list(n, relations, order); + + if(order.size() < n) { // TODO: accept cycles in the poset + throw "corresponding DAG is not acyclic"; + } return relations; } @@ -96,34 +131,8 @@ class Poset { std::vector topologically_sorted_list() const { - std::vector > adjList(n); - std::vector indegree(n, 0); - - for(auto x: order_relations) { - adjList[x.first].push_back(x.second); - indegree[x.second] += 1; - } - - std::queue q; - for(unsigned int i=0; i res; - while(!q.empty()) { - unsigned int curr = q.front(); - res.push_back(curr); - q.pop(); - - for(unsigned int i=0; i + struct parameters { + Point x0; + NT a; + NT eta; + unsigned int order; + NT L; + NT m; + NT kappa; + Eigen::Matrix inv_covariance_matrix; + + parameters(Point x0_, NT a_, NT eta_, Eigen::Matrix inv_covariance_matrix_) + : x0(x0_), a(a_), eta(eta_), order(2), + L(a_), m(a_), kappa(1), + inv_covariance_matrix(inv_covariance_matrix_) {}; + }; + + template + struct GradientFunctor { + typedef typename Point::FT NT; + typedef std::vector pts; + + parameters& params; + + GradientFunctor(parameters& params_) : params(params_) {}; + + Point operator()(unsigned int const& i, pts const& xs, NT const& t) const { + if (i == params.order - 1) { + Point y = xs[0] - params.x0; + Eigen::Matrix result = -2.0 * params.a * params.inv_covariance_matrix * y.getCoefficients(); + return Point(result); + } else { + return xs[i + 1]; + } + } + Point operator()(Point const& x) { + Point y = x - params.x0; + Eigen::Matrix result = -2.0 * params.a * params.inv_covariance_matrix * y.getCoefficients(); + return Point(result); + } + }; + + template + struct FunctionFunctor { + typedef typename Point::FT NT; + + parameters& params; + + FunctionFunctor(parameters& params_) : params(params_) {}; + + NT operator()(Point const& x) const { + Point y = x - params.x0; + Eigen::Matrix y_coeffs = y.getCoefficients(); + return params.a * y_coeffs.dot(params.inv_covariance_matrix * y_coeffs); + } + }; + + template + struct HessianFunctor { + typedef typename Point::FT NT; + + parameters& params; + + HessianFunctor(parameters& params_) : params(params_) {}; + + Point operator()(Point const& x) const { + Eigen::Matrix result = 2.0 * params.a * params.inv_covariance_matrix; + return Point(result); + } + }; +}; + #endif diff --git a/src/volesti/include/preprocess/crhmc/opts.h b/src/volesti/include/preprocess/crhmc/opts.h index 573a69db..fffd29c5 100644 --- a/src/volesti/include/preprocess/crhmc/opts.h +++ b/src/volesti/include/preprocess/crhmc/opts.h @@ -44,6 +44,7 @@ template class opts { bool DynamicWeight = true; //Enable the use of dynamic weights for each variable when sampling bool DynamicStepSize = true; // Enable adaptive step size that avoids low acceptance probability bool DynamicRegularizer = true; //Enable the addition of a regularization term + bool etaInitialize = true; // enable eta initialization Type regularization_factor=1e-20; /*Dynamic step choices*/ Type warmUpStep = 10; diff --git a/src/volesti/include/preprocess/estimate_L_smooth_parameter.hpp b/src/volesti/include/preprocess/estimate_L_smooth_parameter.hpp index 180f2094..340ece75 100644 --- a/src/volesti/include/preprocess/estimate_L_smooth_parameter.hpp +++ b/src/volesti/include/preprocess/estimate_L_smooth_parameter.hpp @@ -21,7 +21,7 @@ template typename NegativeGradientFunctor, typename RandomNumberGenerator > -double estimate_L_smooth(Polytope &P, Point &p, unsigned int const& walk_length, +double estimate_L_smooth(Polytope &P, Point &p, unsigned int const& walk_length, NegativeGradientFunctor F, RandomNumberGenerator &rng) { typedef typename Point::FT NT; @@ -44,11 +44,11 @@ double estimate_L_smooth(Polytope &P, Point &p, unsigned int const& walk_length, RandomWalk walk(P, p, rng); for (unsigned int i=0; i::lowest(), Ltemp; for (auto pit=listOfPoints.begin(); pit!=(listOfPoints.end()-1); ++pit) diff --git a/src/volesti/include/preprocess/inscribed_ellipsoid_rounding.hpp b/src/volesti/include/preprocess/inscribed_ellipsoid_rounding.hpp index dd20d533..4b027f61 100644 --- a/src/volesti/include/preprocess/inscribed_ellipsoid_rounding.hpp +++ b/src/volesti/include/preprocess/inscribed_ellipsoid_rounding.hpp @@ -12,14 +12,9 @@ #define INSCRIBED_ELLIPSOID_ROUNDING_HPP #include "preprocess/max_inscribed_ellipsoid.hpp" -#include "preprocess/analytic_center_linear_ineq.h" +#include "preprocess/barrier_center_ellipsoid.hpp" #include "preprocess/feasible_point.hpp" -enum EllipsoidType -{ - MAX_ELLIPSOID = 1, - LOG_BARRIER = 2 -}; template inline static std::tuple @@ -30,9 +25,11 @@ compute_inscribed_ellipsoid(Custom_MT A, VT b, VT const& x0, if constexpr (ellipsoid_type == EllipsoidType::MAX_ELLIPSOID) { return max_inscribed_ellipsoid(A, b, x0, maxiter, tol, reg); - } else if constexpr (ellipsoid_type == EllipsoidType::LOG_BARRIER) + } else if constexpr (ellipsoid_type == EllipsoidType::LOG_BARRIER || + ellipsoid_type == EllipsoidType::VOLUMETRIC_BARRIER || + ellipsoid_type == EllipsoidType::VAIDYA_BARRIER) { - return analytic_center_linear_ineq(A, b, x0); + return barrier_center_ellipsoid_linear_ineq(A, b, x0); } else { std::runtime_error("Unknown rounding method."); diff --git a/src/volesti/include/preprocess/max_inscribed_ball.hpp b/src/volesti/include/preprocess/max_inscribed_ball.hpp index 573813dd..5f95970b 100644 --- a/src/volesti/include/preprocess/max_inscribed_ball.hpp +++ b/src/volesti/include/preprocess/max_inscribed_ball.hpp @@ -11,7 +11,7 @@ #ifndef MAX_INSCRIBED_BALL_HPP #define MAX_INSCRIBED_BALL_HPP -#include "preprocess/mat_computational_operators.h" +#include "preprocess/rounding_util_functions.hpp" /* This implmentation computes the largest inscribed ball in a given convex polytope P. diff --git a/src/volesti/include/preprocess/max_inscribed_ellipsoid.hpp b/src/volesti/include/preprocess/max_inscribed_ellipsoid.hpp index f44f60ff..cc65006b 100644 --- a/src/volesti/include/preprocess/max_inscribed_ellipsoid.hpp +++ b/src/volesti/include/preprocess/max_inscribed_ellipsoid.hpp @@ -15,7 +15,7 @@ #include #include -#include "preprocess/mat_computational_operators.h" +#include "preprocess/rounding_util_functions.hpp" /* diff --git a/src/volesti/include/random_walks/compute_diameter.hpp b/src/volesti/include/random_walks/compute_diameter.hpp index b20aa839..71007d9e 100644 --- a/src/volesti/include/random_walks/compute_diameter.hpp +++ b/src/volesti/include/random_walks/compute_diameter.hpp @@ -30,11 +30,11 @@ struct compute_diameter }; -template -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/src/volesti/include/random_walks/crhmc/additional_units/dynamic_step_size.hpp b/src/volesti/include/random_walks/crhmc/additional_units/dynamic_step_size.hpp index 36721196..3daee1c3 100644 --- a/src/volesti/include/random_walks/crhmc/additional_units/dynamic_step_size.hpp +++ b/src/volesti/include/random_walks/crhmc/additional_units/dynamic_step_size.hpp @@ -50,7 +50,7 @@ class dynamic_step_size { consecutiveBadStep = IVT::Zero(simdLen); rejectSinceShrink = VT::Zero(simdLen); - if (options.warmUpStep > 0) { + if (options.warmUpStep > 0 && options.etaInitialize) { eta = 1e-3; } else { warmupFinished = true; diff --git a/src/volesti/include/random_walks/crhmc/crhmc_walk.hpp b/src/volesti/include/random_walks/crhmc/crhmc_walk.hpp index 21ed8d6b..9f1dd24a 100644 --- a/src/volesti/include/random_walks/crhmc/crhmc_walk.hpp +++ b/src/volesti/include/random_walks/crhmc/crhmc_walk.hpp @@ -164,6 +164,9 @@ struct CRHMCWalk { inline void disable_adaptive(){ update_modules=false; } + NT get_current_eta() const { + return solver->get_eta(); + } inline void apply(RandomNumberGenerator &rng, int walk_length = 1, bool metropolis_filter = true) diff --git a/src/volesti/include/random_walks/gaussian_accelerated_billiard_walk.hpp b/src/volesti/include/random_walks/gaussian_accelerated_billiard_walk.hpp index fbf43f44..532d7672 100644 --- a/src/volesti/include/random_walks/gaussian_accelerated_billiard_walk.hpp +++ b/src/volesti/include/random_walks/gaussian_accelerated_billiard_walk.hpp @@ -3,14 +3,18 @@ // Copyright (c) 2012-2020 Vissarion Fisikopoulos // Copyright (c) 2018-2020 Apostolos Chalkis // Copyright (c) 2021 Vaibhav Thakkar +// Copyright (c) 2024 Luca Perju // Contributed and/or modified by Vaibhav Thakkar, as part of Google Summer of Code 2021 program. +// Contributed and/or modified by Luca Perju, as part of Google Summer of Code 2024 program. // Licensed under GNU LGPL.3, see LICENCE file #ifndef RANDOM_WALKS_GAUSSIAN_ACCELERATED_BILLIARD_WALK_HPP #define RANDOM_WALKS_GAUSSIAN_ACCELERATED_BILLIARD_WALK_HPP +#include +#include #include "convex_bodies/orderpolytope.h" #include "convex_bodies/ellipsoid.h" #include "convex_bodies/ballintersectconvex.h" @@ -20,8 +24,6 @@ #include "random_walks/compute_diameter.hpp" -// Billiard walk which accelarates each step for uniform distribution and also takes into account -// the shape of the polytope for generating directions. struct GaussianAcceleratedBilliardWalk { @@ -57,56 +59,86 @@ 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; - template + void computeLcov(const E_type E) + { + if constexpr (std::is_base_of, E_type >::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, - Ellipsoid const& E, // ellipsoid representing the Gaussian distribution + E_type const& E, // covariance matrix representing the Gaussian distribution RandomNumberGenerator &rng) { + if(!P.is_normalized()) { + P.normalize(); + } _update_parameters = update_parameters(); - _Len = compute_diameter - ::template compute(P); - - // Removed as will be used for sparse matrices only - // _AA.noalias() = P.get_mat() * P.get_mat().transpose(); - initialize(P, p, E, rng); + _L = compute_diameter::template compute(P); + computeLcov(E); + _E = E; + _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); } - template + template Walk(GenericPolytope& P, Point const& p, - Ellipsoid const& E, // ellipsoid representing the Gaussian distribution + E_type const& E, // covariance matrix representing the Gaussian distribution RandomNumberGenerator &rng, parameters const& params) { + if(!P.is_normalized()) { + P.normalize(); + } _update_parameters = update_parameters(); - _Len = params.set_L ? params.m_L + _L = params.set_L ? params.m_L : compute_diameter ::template compute(P); - - // Removed as will be used for sparse matrices only - // _AA.noalias() = P.get_mat() * P.get_mat().transpose(); - initialize(P, p, E, rng); + computeLcov(E); + _E = E; + _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); } - template - < - typename GenericPolytope, - typename Ellipsoid - > + template inline void apply(GenericPolytope& P, Point &p, // a point to return the result - Ellipsoid const& E, // ellipsoid representing the Gaussian distribution unsigned int const& walk_length, RandomNumberGenerator &rng) { @@ -114,25 +146,50 @@ struct GaussianAcceleratedBilliardWalk NT T; const NT dl = 0.995; int it; + NT coef = 1.0; + NT vEv; for (auto j=0u; j::apply(n, E, rng); - NT norm_v = _v.length(); - _v /= norm_v; + T = -std::log(rng.sample_urdist()) * _L; + + _v = GetDirection::apply(n, rng, false); + _v = Point(_L_cov.template triangularView() * _v.getCoefficients()); + coef = 1.0; + + vEv = (_v.getCoefficients().transpose() * _E.template selfadjointView()).dot(_v.getCoefficients()); Point p0 = _p; - Point v0 = _v; - typename Point::Coeff lambdas0 = _lambdas; - typename Point::Coeff Av0 = _Av; - NT lambda_prev0 = _lambda_prev; it = 0; - while (it < 100*n) + std::pair pbpair; + if(!was_reset) { + pbpair = P.line_positive_intersect(_p, _v, _lambdas, _Av, _lambda_prev, _update_parameters); + } else { + pbpair = P.line_first_positive_intersect(_p, _v, _lambdas, _Av, _update_parameters); + was_reset = false; + } + + if (T <= pbpair.first) { + _p += (T * _v); + _lambda_prev = T; + continue; + } + + _lambda_prev = dl * pbpair.first; + _p += (_lambda_prev * _v); + T -= _lambda_prev; + coef = P.compute_reflection(_v, _p, _AE, _AEA, vEv, _update_parameters); + it++; + + while (it < _rho) { std::pair pbpair - = P.line_positive_intersect(_p, _v, _lambdas, _Av, _lambda_prev, _update_parameters); + = P.line_positive_intersect(_p, _v, _lambdas, _Av, _lambda_prev, _AA, _update_parameters); + _Av *= coef; + _update_parameters.inner_vi_ak *= coef; + pbpair.first /= coef; + if (T <= pbpair.first) { _p += (T * _v); _lambda_prev = T; @@ -141,63 +198,48 @@ struct GaussianAcceleratedBilliardWalk _lambda_prev = dl * pbpair.first; _p += (_lambda_prev * _v); T -= _lambda_prev; - P.compute_reflection(_v, _p, _update_parameters); + coef = P.compute_reflection(_v, _p, _AE, _AEA, vEv, _update_parameters); it++; } - if (it == 100*n) - { + if (it == _rho){ _p = p0; - _lambdas = lambdas0; - _Av = Av0; - _lambda_prev = lambda_prev0; - } - else - { - // metropolis filter - NT u_logprob = log(rng.sample_urdist()); - NT accept_logprob = 0.5 * (norm_v * norm_v) * (E.mat_mult(v0) - E.mat_mult(_v)); - // std::cout << "diff: " << (accept_logprob - u_logprob) << std::endl; - if(u_logprob > accept_logprob) { - // reject - _p = p0; - _lambdas = lambdas0; - _Av = Av0; - _lambda_prev = lambda_prev0; - } + was_reset = true; } } - p = _p; - } - inline void update_delta(NT L) - { - _Len = L; + p = _p; } private : - template - < - typename GenericPolytope, - typename Ellipsoid - > + template inline void initialize(GenericPolytope& P, Point const& p, // a point to start - Ellipsoid const& E, // ellipsoid representing the Gaussian distribution RandomNumberGenerator &rng) { unsigned int n = P.dimension(); const NT dl = 0.995; + was_reset = false; _lambdas.setZero(P.num_of_hyperplanes()); _Av.setZero(P.num_of_hyperplanes()); _p = p; - _v = GetGaussianDirection::apply(n, E, rng); - NT norm_v = _v.length(); - _v /= norm_v; + _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)); + }*/ - NT T = -std::log(rng.sample_urdist()) * _Len; + _v = GetDirection::apply(n, rng, false); + _v = Point(_L_cov.template triangularView() * _v.getCoefficients()); + + NT T = -std::log(rng.sample_urdist()) * _L; Point p0 = _p; int it = 0; + NT coef = 1.0; + NT vEv = (_v.getCoefficients().transpose() * _E.template selfadjointView()).dot(_v.getCoefficients()); std::pair pbpair = P.line_first_positive_intersect(_p, _v, _lambdas, _Av, _update_parameters); @@ -209,17 +251,21 @@ struct GaussianAcceleratedBilliardWalk _lambda_prev = dl * pbpair.first; _p += (_lambda_prev * _v); T -= _lambda_prev; - P.compute_reflection(_v, _p, _update_parameters); + coef = P.compute_reflection(_v, _p, _AE, _AEA, vEv, _update_parameters); - while (it <= 100*n) + while (it <= _rho) { std::pair pbpair - = P.line_positive_intersect(_p, _v, _lambdas, _Av, _lambda_prev, _update_parameters); + = P.line_positive_intersect(_p, _v, _lambdas, _Av, _lambda_prev, _AA, _update_parameters); + _Av *= coef; + _update_parameters.inner_vi_ak *= coef; + pbpair.first /= coef; + if (T <= pbpair.first) { _p += (T * _v); _lambda_prev = T; break; - } else if (it == 100*n) { + } else if (it == _rho) { _lambda_prev = rng.sample_urdist() * pbpair.first; _p += (_lambda_prev * _v); break; @@ -227,19 +273,25 @@ struct GaussianAcceleratedBilliardWalk _lambda_prev = dl * pbpair.first; _p += (_lambda_prev * _v); T -= _lambda_prev; - P.compute_reflection(_v, _p, _update_parameters); + coef = P.compute_reflection(_v, _p, _AE, _AEA, vEv, _update_parameters); it++; } } - NT _Len; + NT _L; Point _p; Point _v; NT _lambda_prev; - // MT _AA; // Removed as will be used for sparse matrices only + DenseMT _AA; + E_type _L_cov; // LL' = inv(E) + DenseMT _AE; + E_type _E; + VT _AEA; + unsigned int _rho; update_parameters _update_parameters; typename Point::Coeff _lambdas; typename Point::Coeff _Av; + bool was_reset; }; }; diff --git a/src/volesti/include/random_walks/gaussian_helpers.hpp b/src/volesti/include/random_walks/gaussian_helpers.hpp index d1ab5429..76a15113 100644 --- a/src/volesti/include/random_walks/gaussian_helpers.hpp +++ b/src/volesti/include/random_walks/gaussian_helpers.hpp @@ -9,7 +9,14 @@ NT eval_exp(Point const& p, NT const& a) { return std::exp(-a * p.squared_length()); } - +template +NT eval_exp(Point const& p, MT const& inv_covariance_matrix, NT const& a_next, NT const& a_curr) +{ + Eigen::Matrix dist_vector = p.getCoefficients(); + NT mahalanobis_dist = dist_vector.transpose() * inv_covariance_matrix * dist_vector; + NT log_ratio = (a_curr - a_next) * mahalanobis_dist; + return std::exp(log_ratio); +} template NT get_max(Point const& l, Point const& u, NT const& a_i) diff --git a/src/volesti/include/random_walks/uniform_accelerated_billiard_walk.hpp b/src/volesti/include/random_walks/uniform_accelerated_billiard_walk.hpp index 980c8e5d..f91abee4 100644 --- a/src/volesti/include/random_walks/uniform_accelerated_billiard_walk.hpp +++ b/src/volesti/include/random_walks/uniform_accelerated_billiard_walk.hpp @@ -4,6 +4,7 @@ // Copyright (c) 2018-2020 Apostolos Chalkis // Contributed and/or modified by Alexandros Manochis, as part of Google Summer of Code 2020 program. +// Contributed and/or modified by Luca Perju, as part of Google Summer of Code 2024 program. // Licensed under GNU LGPL.3, see LICENCE file @@ -11,6 +12,126 @@ #define RANDOM_WALKS_ACCELERATED_IMPROVED_BILLIARD_WALK_HPP #include "sampling/sphere.hpp" +#include +#include +#include + +const double eps = 1e-10; + +// data structure which maintains the values of (b - Ar)/Av, and can extract the minimum positive value and the facet associated with it +// vec[i].first contains the value of (b(i) - Ar(i))/Av(i) + moved_dist, where moved_dist is the total distance that the point has travelled so far +// The heap will only contain the values from vec which are greater than moved_dist (so that they are positive) +template +class BoundaryOracleHeap { +public: + int n, heap_size; + std::vector> heap; + std::vector> vec; + +private: + int siftDown(int index) { + while((index << 1) + 1 < heap_size) { + int child = (index << 1) + 1; + if(child + 1 < heap_size && heap[child + 1].first < heap[child].first - eps) { + child += 1; + } + if(heap[child].first < heap[index].first - eps) + { + std::swap(heap[child], heap[index]); + std::swap(vec[heap[child].second].second, vec[heap[index].second].second); + index = child; + } else { + return index; + } + } + return index; + } + + int siftUp(int index) { + while(index > 0 && heap[(index - 1) >> 1].first - eps > heap[index].first) { + std::swap(heap[(index - 1) >> 1], heap[index]); + std::swap(vec[heap[(index - 1) >> 1].second].second, vec[heap[index].second].second); + index = (index - 1) >> 1; + } + return index; + } + + // takes the index of a facet, and (in case it is in the heap) removes it from the heap. + void remove (int index) { + index = vec[index].second; + if(index == -1) { + return; + } + std::swap(heap[heap_size - 1], heap[index]); + std::swap(vec[heap[heap_size - 1].second].second, vec[heap[index].second].second); + vec[heap[heap_size - 1].second].second = -1; + heap_size -= 1; + index = siftDown(index); + siftUp(index); + } + + // inserts a new value into the heap, with its associated facet + void insert (const std::pair val) { + vec[val.second].second = heap_size; + vec[val.second].first = val.first; + heap[heap_size++] = val; + siftUp(heap_size - 1); + } + +public: + BoundaryOracleHeap() {} + + BoundaryOracleHeap(int n) : n(n), heap_size(0) { + heap.resize(n); + vec.resize(n); + } + + // rebuilds the heap with the existing values from vec + // O(n) + void rebuild (const NT &moved_dist) { + heap_size = 0; + for(int i = 0; i < n; ++i) { + vec[i].second = -1; + if(vec[i].first - eps > moved_dist) { + vec[i].second = heap_size; + heap[heap_size++] = {vec[i].first, i}; + } + } + for(int i = heap_size - 1; i >= 0; --i) { + siftDown(i); + } + } + + // returns (b(i) - Ar(i))/Av(i) + moved_dist + // O(1) + NT get_val (const int &index) { + return vec[index].first; + } + + // returns the nearest facet + // O(1) + std::pair get_min () { + return heap[0]; + } + + // changes the stored value for a given facet, and updates the heap accordingly + // O(logn) + void change_val(const int& index, const NT& new_val, const NT& moved_dist) { + if(new_val < moved_dist - eps) { + vec[index].first = new_val; + remove(index); + } else { + if(vec[index].second == -1) { + insert({new_val, index}); + } else { + heap[vec[index].second].first = new_val; + vec[index].first = new_val; + siftDown(vec[index].second); + siftUp(vec[index].second); + } + } + } +}; // Billiard walk which accelarates each step for uniform distribution @@ -37,35 +158,46 @@ struct AcceleratedBilliardWalk struct update_parameters { update_parameters() - : facet_prev(0), hit_ball(false), inner_vi_ak(0.0), ball_inner_norm(0.0) + : facet_prev(0), hit_ball(false), inner_vi_ak(0.0), ball_inner_norm(0.0), moved_dist(0.0) {} int facet_prev; bool hit_ball; double inner_vi_ak; double ball_inner_norm; + double moved_dist; }; parameters param; template - < - typename Polytope, - typename RandomNumberGenerator - > + < + typename Polytope, + typename RandomNumberGenerator + > struct Walk { typedef typename Polytope::PointType Point; typedef typename Polytope::MT MT; + typedef typename Eigen::Matrix DenseMT; typedef typename Point::FT NT; + using AA_type = std::conditional_t< std::is_same_v>, typename Eigen::SparseMatrix, DenseMT >; + // AA is sparse colMajor if MT is sparse rowMajor, and Dense otherwise template Walk(GenericPolytope &P, Point const& p, RandomNumberGenerator &rng) { + if(!P.is_normalized()) { + P.normalize(); + } _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() = (DenseMT)(P.get_mat() * P.get_mat().transpose()); + } _rho = 1000 * P.dimension(); // upper bound for the number of reflections (experimental) initialize(P, p, rng); } @@ -74,11 +206,18 @@ struct AcceleratedBilliardWalk Walk(GenericPolytope &P, Point const& p, RandomNumberGenerator &rng, parameters const& params) { + if(!P.is_normalized()) { + P.normalize(); + } _update_parameters = update_parameters(); _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() = (DenseMT)(P.get_mat() * P.get_mat().transpose()); + } _rho = 1000 * P.dimension(); // upper bound for the number of reflections (experimental) initialize(P, p, rng); } @@ -96,6 +235,12 @@ struct AcceleratedBilliardWalk NT T; const NT dl = 0.995; int it; + typename Point::Coeff b; + NT* b_data; + if constexpr (std::is_same>::value) { + b = P.get_vec(); + b_data = b.data(); + } for (auto j=0u; j pbpair = P.line_positive_intersect(_p, _v, _lambdas, _Av, - _lambda_prev, _update_parameters); + std::pair pbpair = P.line_first_positive_intersect(_p, _v, _lambdas, _Av, _update_parameters); + if (T <= pbpair.first) { _p += (T * _v); _lambda_prev = T; @@ -113,28 +258,52 @@ struct AcceleratedBilliardWalk } _lambda_prev = dl * pbpair.first; - _p += (_lambda_prev * _v); + if constexpr (std::is_same>::value) { + _update_parameters.moved_dist = _lambda_prev; + NT* Ar_data = _lambdas.data(); + NT* Av_data = _Av.data(); + for(int i = 0; i < P.num_of_hyperplanes(); ++i) { + _distances_set.vec[i].first = ( *(b_data + i) - (*(Ar_data + i)) ) / (*(Av_data + i)); + } + // rebuild the heap with the new values of (b - Ar) / Av + _distances_set.rebuild(_update_parameters.moved_dist); + } else { + _p += (_lambda_prev * _v); + } T -= _lambda_prev; P.compute_reflection(_v, _p, _update_parameters); it++; while (it < _rho) - { - std::pair pbpair - = P.line_positive_intersect(_p, _v, _lambdas, _Av, _lambda_prev, - _AA, _update_parameters); + { + std::pair pbpair; + if constexpr (std::is_same>::value) { + pbpair = P.line_positive_intersect(_p, _lambdas, _Av, _lambda_prev, + _distances_set, _AA, _update_parameters); + } else { + pbpair = P.line_positive_intersect(_p, _v, _lambdas, _Av, _lambda_prev, + _AA, _update_parameters); + } if (T <= pbpair.first) { _p += (T * _v); _lambda_prev = T; break; } _lambda_prev = dl * pbpair.first; - _p += (_lambda_prev * _v); + if constexpr (std::is_same>::value) { + _update_parameters.moved_dist += _lambda_prev; + } else { + _p += (_lambda_prev * _v); + } T -= _lambda_prev; P.compute_reflection(_v, _p, _update_parameters); it++; } - if (it == _rho) _p = p0; + _p += _update_parameters.moved_dist * _v; + _update_parameters.moved_dist = 0.0; + if (it == _rho) { + _p = p0; + } } p = _p; } @@ -185,8 +354,7 @@ struct AcceleratedBilliardWalk apply(P, p, walk_length, rng); max_dist = get_max_distance(pointset, p, rad); - if (max_dist > Lmax) - { + if (max_dist > Lmax) { Lmax = max_dist; } if (2.0*rad > Lmax) { @@ -196,8 +364,7 @@ struct AcceleratedBilliardWalk } if (Lmax > _L) { - if (P.dimension() <= 2500) - { + if (P.dimension() <= 2500) { update_delta(Lmax); } else{ @@ -235,6 +402,7 @@ struct AcceleratedBilliardWalk _Av.setZero(P.num_of_hyperplanes()); _p = p; _v = GetDirection::apply(n, rng); + _distances_set = BoundaryOracleHeap(P.num_of_hyperplanes()); NT T = -std::log(rng.sample_urdist()) * _L; Point p0 = _p; @@ -296,11 +464,12 @@ struct AcceleratedBilliardWalk Point _p; Point _v; NT _lambda_prev; - MT _AA; + AA_type _AA; unsigned int _rho; update_parameters _update_parameters; typename Point::Coeff _lambdas; typename Point::Coeff _Av; + BoundaryOracleHeap _distances_set; }; }; diff --git a/src/volesti/include/random_walks/uniform_accelerated_billiard_walk_parallel.hpp b/src/volesti/include/random_walks/uniform_accelerated_billiard_walk_parallel.hpp index 1c2630f9..7ec437b2 100644 --- a/src/volesti/include/random_walks/uniform_accelerated_billiard_walk_parallel.hpp +++ b/src/volesti/include/random_walks/uniform_accelerated_billiard_walk_parallel.hpp @@ -53,10 +53,10 @@ struct AcceleratedBilliardWalkParallel template - < - typename Polytope, - typename RandomNumberGenerator - > + < + typename Polytope, + typename RandomNumberGenerator + > struct Walk { typedef typename Polytope::PointType Point; @@ -66,6 +66,9 @@ struct AcceleratedBilliardWalkParallel template Walk(GenericPolytope &P) { + if(!P.is_normalized()) { + P.normalize(); + } _L = compute_diameter ::template compute(P); _AA.noalias() = P.get_mat() * P.get_mat().transpose(); @@ -76,6 +79,9 @@ struct AcceleratedBilliardWalkParallel template Walk(GenericPolytope &P, NT const& L) { + if(!P.is_normalized()) { + P.normalize(); + } _L = L > NT(0) ? L : compute_diameter ::template compute(P); @@ -190,8 +196,7 @@ struct AcceleratedBilliardWalkParallel apply(P, params, walk_length, rng); max_dist = get_max_distance(pointset, params.p, rad); - if (max_dist > Lmax) - { + if (max_dist > Lmax) { Lmax = max_dist; } if (2.0*rad > Lmax) { @@ -200,10 +205,8 @@ struct AcceleratedBilliardWalkParallel pointset.push_back(params.p); } - if (Lmax > _L) - { - if (P.dimension() <= 2500) - { + if (Lmax > _L) { + if (P.dimension() <= 2500) { update_delta(Lmax); } else{ diff --git a/src/volesti/include/sampling/mmcs.hpp b/src/volesti/include/sampling/mmcs.hpp index 07529ad6..02791535 100644 --- a/src/volesti/include/sampling/mmcs.hpp +++ b/src/volesti/include/sampling/mmcs.hpp @@ -81,7 +81,7 @@ bool perform_mmcs_step(Polytope &P, walk.template get_starting_point(P, p, q, 10, rng); for (int i = 0; i < window; i++) { - walk.template apply(P, q, walk_length, rng); + walk.apply(P, q, walk_length, rng); winPoints.col(i) = q.getCoefficients(); } estimator.update_estimator(winPoints); diff --git a/src/volesti/include/sampling/parallel_mmcs.hpp b/src/volesti/include/sampling/parallel_mmcs.hpp index a8308f4c..da22e8f6 100644 --- a/src/volesti/include/sampling/parallel_mmcs.hpp +++ b/src/volesti/include/sampling/parallel_mmcs.hpp @@ -126,7 +126,7 @@ bool perform_parallel_mmcs_step(Polytope &P, walk.template get_starting_point(P, p, thread_random_walk_parameters, 10, rng); for (int i = 0; i < window; i++) { - walk.template apply(P, thread_random_walk_parameters, walk_length, rng); + walk.apply(P, thread_random_walk_parameters, walk_length, rng); winPoints_per_thread[thread_index].col(i) = thread_random_walk_parameters.p.getCoefficients(); } diff --git a/src/volesti/include/sampling/random_point_generators.hpp b/src/volesti/include/sampling/random_point_generators.hpp index e70e1461..998d8daf 100644 --- a/src/volesti/include/sampling/random_point_generators.hpp +++ b/src/volesti/include/sampling/random_point_generators.hpp @@ -35,7 +35,7 @@ struct RandomPointGenerator Walk walk(P, p, rng, parameters); for (unsigned int i=0; i struct policy_storing { template - static void store(WalkPolicy &policy, PointList &randPoints, ThreadParameters &thread_random_walk_parameters) + static void store(WalkPolicy &policy, PointList &randPoints, ThreadParameters &thread_random_walk_parameters) { policy.apply(randPoints, thread_random_walk_parameters.p); } @@ -30,7 +30,7 @@ template <> struct policy_storing { template - static void store(WalkPolicy &policy, PointList &randPoints, ThreadParameters &thread_random_walk_parameters) + static void store(WalkPolicy &policy, PointList &randPoints, ThreadParameters &thread_random_walk_parameters) { policy.apply(randPoints, thread_random_walk_parameters.p1); policy.apply(randPoints, thread_random_walk_parameters.p2); @@ -38,7 +38,7 @@ struct policy_storing }; template <> -struct policy_storing : policy_storing +struct policy_storing : policy_storing {}; @@ -72,7 +72,7 @@ struct RandomPointGeneratorMultiThread omp_set_num_threads(num_threads); unsigned int d = P.dimension(), m = P.num_of_hyperplanes(); - + std::vector num_points_per_thread(rnum%num_threads, rnum/num_threads+1); std::vector a(num_threads - rnum%num_threads, rnum/num_threads); num_points_per_thread.insert(num_points_per_thread.end(), a.begin(), a.end()); @@ -88,7 +88,7 @@ struct RandomPointGeneratorMultiThread for (unsigned int it = 0; it < num_points_per_thread[thread_index]; it++) { - walk.template apply(P, thread_random_walk_parameters, walk_length, rng); + walk.apply(P, thread_random_walk_parameters, walk_length, rng); #pragma omp critical { policy_storing::template store(policy, randPoints, thread_random_walk_parameters); @@ -135,7 +135,7 @@ struct RandomPointGeneratorMultiThread for (unsigned int it = 0; it < num_points_per_thread[thread_index]; it++) { - walk.template apply(P, thread_random_walk_parameters, walk_length, rng); + walk.apply(P, thread_random_walk_parameters, walk_length, rng); #pragma omp critical { policy_storing::template store(policy, randPoints, thread_random_walk_parameters); @@ -178,7 +178,7 @@ struct GaussianPointGeneratorMultiThread omp_set_num_threads(num_threads); unsigned int d = P.dimension(), m = P.num_of_hyperplanes(); - + std::vector num_points_per_thread(rnum%num_threads, rnum/num_threads+1); std::vector a(num_threads - rnum%num_threads, rnum/num_threads); num_points_per_thread.insert(num_points_per_thread.end(), a.begin(), a.end()); @@ -194,7 +194,7 @@ struct GaussianPointGeneratorMultiThread for (unsigned int it = 0; it < num_points_per_thread[thread_index]; it++) { - walk.template apply(P, thread_random_walk_parameters, a_i, walk_length, rng); + walk.apply(P, thread_random_walk_parameters, a_i, walk_length, rng); #pragma omp critical { policy_storing::template store(policy, randPoints, thread_random_walk_parameters); @@ -226,7 +226,7 @@ struct GaussianPointGeneratorMultiThread omp_set_num_threads(num_threads); unsigned int d = P.dimension(), m = P.num_of_hyperplanes(); - + std::vector num_points_per_thread(rnum%num_threads, rnum/num_threads+1); std::vector a(num_threads - rnum%num_threads, rnum/num_threads); num_points_per_thread.insert(num_points_per_thread.end(), a.begin(), a.end()); @@ -242,7 +242,7 @@ struct GaussianPointGeneratorMultiThread for (unsigned int it = 0; it < num_points_per_thread[thread_index]; it++) { - walk.template apply(P, thread_random_walk_parameters, a_i, walk_length, rng); + walk.apply(P, thread_random_walk_parameters, a_i, walk_length, rng); #pragma omp critical { policy_storing::template store(policy, randPoints, thread_random_walk_parameters); diff --git a/src/volesti/include/sampling/sample_correlation_matrices.hpp b/src/volesti/include/sampling/sample_correlation_matrices.hpp index 50b7a45d..b2bd118b 100644 --- a/src/volesti/include/sampling/sample_correlation_matrices.hpp +++ b/src/volesti/include/sampling/sample_correlation_matrices.hpp @@ -14,6 +14,21 @@ #include +template +MT getCoefficientsFromMatrix(const MT& mat) { + int n = mat.rows(); + int d = n * (n - 1) / 2; + Eigen::Matrix coeffs(d); + int k = 0; + for (int i = 0; i < n; ++i) { + for (int j = 0; j < i; ++j) { + coeffs(k) = mat(i, j); + ++k; + } + } + return coeffs; +} + // New implementations for sampling correlation matrices with CorrelationSpectrahedron and CorrelationSpectrahedron_MT template diff --git a/src/volesti/include/volume/volume_cooling_balls.hpp b/src/volesti/include/volume/volume_cooling_balls.hpp index d90604fa..7f3d00d5 100644 --- a/src/volesti/include/volume/volume_cooling_balls.hpp +++ b/src/volesti/include/volume/volume_cooling_balls.hpp @@ -466,7 +466,7 @@ NT estimate_ratio(PolyBall1 &Pb1, do { - walk.template apply(Pb1, p, walk_length, rng); + walk.apply(Pb1, p, walk_length, rng); } while(!estimate_ratio_generic(Pb2, p, error, ratio_parameters)); return NT(ratio_parameters.count_in) / NT(ratio_parameters.tot_count); @@ -674,14 +674,14 @@ NT estimate_ratio_interval(PolyBall1 &Pb1, for (int i = 0; i < ratio_parameters.W; ++i) { - walk.template apply(Pb1, p, walk_length, rng); + walk.apply(Pb1, p, walk_length, rng); full_sliding_window(Pb2, p, ratio_parameters); } ratio_parameters.mean = ratio_parameters.sum / NT(ratio_parameters.W); do { - walk.template apply(Pb1, p, walk_length, rng); + walk.apply(Pb1, p, walk_length, rng); } while (!estimate_ratio_interval_generic(Pb2, p, error, zp, ratio_parameters)); return NT(ratio_parameters.count_in) / NT(ratio_parameters.tot_count); diff --git a/src/volesti/include/volume/volume_cooling_gaussians.hpp b/src/volesti/include/volume/volume_cooling_gaussians.hpp index ca06ce31..62835bd8 100644 --- a/src/volesti/include/volume/volume_cooling_gaussians.hpp +++ b/src/volesti/include/volume/volume_cooling_gaussians.hpp @@ -241,7 +241,7 @@ void compute_annealing_schedule(Polytope& P, // Compute some ratios to decide if this is the last gaussian for (unsigned int j = 0; j < steps; j++) { - walk.template apply(P, p, a_vals[it], walk_length, rng); + walk.apply(P, p, a_vals[it], walk_length, rng); curr_its += 1.0; curr_fn += eval_exp(p, next_a) / eval_exp(p, a_vals[it]); } @@ -398,7 +398,7 @@ double volume_cooling_gaussians(Polytope& Pin, while (!done || (*itsIt)