From 6d73fb516ea02fbebd5b5dacda23e7e106fdff04 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 +++ .../preprocess/analytic_center_linear_ineq.h | 134 ---- .../preprocess/barrier_center_ellipsoid.hpp | 120 ++++ 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 +- ...perators.h => rounding_util_functions.hpp} | 146 ++++- .../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 +- .../volume/volume_cooling_gaussians_crhmc.hpp | 459 ++++++++++++++ ...e_cooling_nonspherical_gaussians_crhmc.hpp | 594 ++++++++++++++++++ 40 files changed, 2123 insertions(+), 392 deletions(-) delete mode 100644 src/volesti/include/preprocess/analytic_center_linear_ineq.h create mode 100644 src/volesti/include/preprocess/barrier_center_ellipsoid.hpp rename src/volesti/include/preprocess/{mat_computational_operators.h => rounding_util_functions.hpp} (65%) create mode 100644 src/volesti/include/volume/volume_cooling_gaussians_crhmc.hpp create mode 100644 src/volesti/include/volume/volume_cooling_nonspherical_gaussians_crhmc.hpp 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/analytic_center_linear_ineq.h b/src/volesti/include/preprocess/analytic_center_linear_ineq.h deleted file mode 100644 index 80d731f6..00000000 --- a/src/volesti/include/preprocess/analytic_center_linear_ineq.h +++ /dev/null @@ -1,134 +0,0 @@ -// VolEsti (volume computation and sampling library) - -// Copyright (c) 2024 Vissarion Fisikopoulos -// Copyright (c) 2024 Apostolos Chalkis -// Copyright (c) 2024 Elias Tsigaridas - -// Licensed under GNU LGPL.3, see LICENCE file - - -#ifndef ANALYTIC_CENTER_H -#define ANALYTIC_CENTER_H - -#include - -#include "preprocess/max_inscribed_ball.hpp" -#include "preprocess/feasible_point.hpp" -#include "preprocess/mat_computational_operators.h" - -template -NT get_max_step(VT const& Ad, VT const& b_Ax) -{ - const int m = Ad.size(); - NT max_element = std::numeric_limits::lowest(), max_element_temp; - for (int i = 0; i < m; i++) { - max_element_temp = Ad.coeff(i) / b_Ax.coeff(i); - if (max_element_temp > max_element) { - max_element = max_element_temp; - } - } - - return NT(1) / max_element; -} - -template -void get_hessian_grad_logbarrier(MT const& A, MT const& A_trans, VT const& b, - VT const& x, VT const& Ax, MT &H, VT &grad, VT &b_Ax) -{ - b_Ax.noalias() = b - Ax; - VT s = b_Ax.cwiseInverse(); - VT s_sq = s.cwiseProduct(s); - // Gradient of the log-barrier function - grad.noalias() = A_trans * s; - // Hessian of the log-barrier function - update_Atrans_Diag_A(H, A_trans, A, s_sq.asDiagonal()); -} - -/* - This implementation computes the analytic center of a polytope given - as a set of linear inequalities P = {x | Ax <= b}. The analytic center - is the minimizer of the log barrier function i.e., the optimal solution - of the following optimization problem (Convex Optimization, Boyd and Vandenberghe, Section 8.5.3), - - \min -\sum \log(b_i - a_i^Tx), where a_i is the i-th row of A. - - The function solves the problem by using the Newton method. - - Input: (i) Matrix A, vector b such that the polytope P = {x | Ax<=b} - (ii) The number of maximum iterations, max_iters - (iii) Tolerance parameter grad_err_tol to bound the L2-norm of the gradient - (iv) Tolerance parameter rel_pos_err_tol to check the relative progress in each iteration - - Output: (i) The Hessian computed on the analytic center - (ii) The analytic center of the polytope - (iii) A boolean variable that declares convergence - - Note: Using MT as to deal with both dense and sparse matrices, MT_dense will be the type of result matrix -*/ -template -std::tuple analytic_center_linear_ineq(MT const& A, VT const& b, VT const& x0, - unsigned int const max_iters = 500, - NT const grad_err_tol = 1e-08, - NT const rel_pos_err_tol = 1e-12) -{ - // Initialization - VT x = x0; - VT Ax = A * x; - const int n = A.cols(), m = A.rows(); - MT H(n, n), A_trans = A.transpose(); - VT grad(n), d(n), Ad(m), b_Ax(m), step_d(n), x_prev; - NT grad_err, rel_pos_err, rel_pos_err_temp, step; - unsigned int iter = 0; - bool converged = false; - const NT tol_bnd = NT(0.01); - - auto llt = initialize_chol(A_trans, A); - get_hessian_grad_logbarrier(A, A_trans, b, x, Ax, H, grad, b_Ax); - - do { - iter++; - // Compute the direction - d.noalias() = - solve_vec(llt, H, grad); - Ad.noalias() = A * d; - // Compute the step length - step = std::min((NT(1) - tol_bnd) * get_max_step(Ad, b_Ax), NT(1)); - step_d.noalias() = step*d; - x_prev = x; - x += step_d; - Ax.noalias() += step*Ad; - - // Compute the max_i\{ |step*d_i| ./ |x_i| \} - rel_pos_err = std::numeric_limits::lowest(); - for (int i = 0; i < n; i++) - { - rel_pos_err_temp = std::abs(step_d.coeff(i) / x_prev.coeff(i)); - if (rel_pos_err_temp > rel_pos_err) - { - rel_pos_err = rel_pos_err_temp; - } - } - - get_hessian_grad_logbarrier(A, A_trans, b, x, Ax, H, grad, b_Ax); - grad_err = grad.norm(); - - if (iter >= max_iters || grad_err <= grad_err_tol || rel_pos_err <= rel_pos_err_tol) - { - converged = true; - break; - } - } while (true); - - return std::make_tuple(MT_dense(H), x, converged); -} - -template -std::tuple analytic_center_linear_ineq(MT const& A, VT const& b, - unsigned int const max_iters = 500, - NT const grad_err_tol = 1e-08, - NT const rel_pos_err_tol = 1e-12) -{ - VT x0 = compute_feasible_point(A, b); - return analytic_center_linear_ineq(A, b, x0, max_iters, grad_err_tol, rel_pos_err_tol); -} - -#endif diff --git a/src/volesti/include/preprocess/barrier_center_ellipsoid.hpp b/src/volesti/include/preprocess/barrier_center_ellipsoid.hpp new file mode 100644 index 00000000..9f440f22 --- /dev/null +++ b/src/volesti/include/preprocess/barrier_center_ellipsoid.hpp @@ -0,0 +1,120 @@ +// VolEsti (volume computation and sampling library) + +// Copyright (c) 2024 Vissarion Fisikopoulos +// Copyright (c) 2024 Apostolos Chalkis +// Copyright (c) 2024 Elias Tsigaridas + +// Licensed under GNU LGPL.3, see LICENCE file + + +#ifndef BARRIER_CENTER_ELLIPSOID_HPP +#define BARRIER_CENTER_ELLIPSOID_HPP + +#include + +#include "preprocess/max_inscribed_ball.hpp" +#include "preprocess/feasible_point.hpp" +#include "preprocess/rounding_util_functions.hpp" + +/* + This implementation computes the analytic or the volumetric center of a polytope given + as a set of linear inequalities P = {x | Ax <= b}. The analytic center is the tminimizer + of the log barrier function, i.e., the optimal solution + of the following optimization problem (Convex Optimization, Boyd and Vandenberghe, Section 8.5.3), + + \min -\sum \log(b_i - a_i^Tx), where a_i is the i-th row of A. + + The volumetric center is the minimizer of the volumetric barrier function, i.e., the optimal + solution of the following optimization problem, + + \min logdet \nabla^2 f(x), where f(x) the log barrier function + + The Vaidya center is the minimizer of the Vaidya barrier function, i.e., the optimal + solution of the following optimization problem, + + \min logdet \nabla^2 f(x) + d/m f(x), where f(x) the log barrier function. + + The function solves the problems by using the Newton method. + + Input: (i) Matrix A, vector b such that the polytope P = {x | Ax<=b} + (ii) The number of maximum iterations, max_iters + (iii) Tolerance parameter grad_err_tol to bound the L2-norm of the gradient + (iv) Tolerance parameter rel_pos_err_tol to check the relative progress in each iteration + + Output: (i) The Hessian of the barrier function + (ii) The analytic/volumetric center of the polytope + (iii) A boolean variable that declares convergence + + Note: Using MT as to deal with both dense and sparse matrices, MT_dense will be the type of result matrix +*/ +template +std::tuple barrier_center_ellipsoid_linear_ineq(MT const& A, VT const& b, VT const& x0, + unsigned int const max_iters = 500, + NT const grad_err_tol = 1e-08, + NT const rel_pos_err_tol = 1e-12) +{ + // Initialization + VT x = x0; + VT Ax = A * x; + const int n = A.cols(), m = A.rows(); + MT H(n, n), A_trans = A.transpose(); + VT grad(n), d(n), Ad(m), b_Ax(m), step_d(n), x_prev; + NT grad_err, rel_pos_err, rel_pos_err_temp, step, obj_val, obj_val_prev; + unsigned int iter = 0; + bool converged = false; + const NT tol_bnd = NT(0.01), tol_obj = NT(1e-06); + + auto [step_iter, max_step_multiplier] = init_step(); + auto llt = initialize_chol(A_trans, A); + get_barrier_hessian_grad(A, A_trans, b, x, Ax, llt, + H, grad, b_Ax, obj_val); + do { + iter++; + // Compute the direction + d.noalias() = - solve_vec(llt, H, grad); + Ad.noalias() = A * d; + // Compute the step length + step = std::min(max_step_multiplier * get_max_step(Ad, b_Ax), step_iter); + step_d.noalias() = step*d; + x_prev = x; + x += step_d; + Ax.noalias() += step*Ad; + + // Compute the max_i\{ |step*d_i| ./ |x_i| \} + rel_pos_err = std::numeric_limits::lowest(); + for (int i = 0; i < n; i++) + { + rel_pos_err_temp = std::abs(step_d.coeff(i) / x_prev.coeff(i)); + if (rel_pos_err_temp > rel_pos_err) + { + rel_pos_err = rel_pos_err_temp; + } + } + + obj_val_prev = obj_val; + get_barrier_hessian_grad(A, A_trans, b, x, Ax, llt, + H, grad, b_Ax, obj_val); + grad_err = grad.norm(); + + if (iter >= max_iters || grad_err <= grad_err_tol || rel_pos_err <= rel_pos_err_tol) + { + converged = true; + break; + } + get_step_next_iteration(obj_val_prev, obj_val, tol_obj, step_iter); + } while (true); + + return std::make_tuple(MT_dense(H), x, converged); +} + +template +std::tuple barrier_center_ellipsoid_linear_ineq(MT const& A, VT const& b, + unsigned int const max_iters = 500, + NT const grad_err_tol = 1e-08, + NT const rel_pos_err_tol = 1e-12) +{ + VT x0 = compute_feasible_point(A, b); + return barrier_center_ellipsoid_linear_ineq(A, b, x0, max_iters, grad_err_tol, rel_pos_err_tol); +} + +#endif // BARRIER_CENTER_ELLIPSOID_HPP 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/preprocess/mat_computational_operators.h b/src/volesti/include/preprocess/rounding_util_functions.hpp similarity index 65% rename from src/volesti/include/preprocess/mat_computational_operators.h rename to src/volesti/include/preprocess/rounding_util_functions.hpp index af200d22..bb07d5b3 100644 --- a/src/volesti/include/preprocess/mat_computational_operators.h +++ b/src/volesti/include/preprocess/rounding_util_functions.hpp @@ -1,14 +1,15 @@ // VolEsti (volume computation and sampling library) -// Copyright (c) 2024 Vissarion Fisikopoulos -// Copyright (c) 2024 Apostolos Chalkis -// Copyright (c) 2024 Elias Tsigaridas +// Copyright (c) 2012-2024 Vissarion Fisikopoulos +// Copyright (c) 2018-2024 Apostolos Chalkis + +//Contributed and/or modified by Alexandros Manochis, as part of Google Summer of Code 2020 program. // Licensed under GNU LGPL.3, see LICENCE file -#ifndef MAT_COMPUTATIONAL_OPERATORS_H -#define MAT_COMPUTATIONAL_OPERATORS_H +#ifndef ROUNDING_UTIL_FUNCTIONS_HPP +#define ROUNDING_UTIL_FUNCTIONS_HPP #include @@ -17,9 +18,36 @@ #include "Spectra/include/Spectra/MatOp/SparseSymMatProd.h" +enum EllipsoidType +{ + MAX_ELLIPSOID = 1, + LOG_BARRIER = 2, + VOLUMETRIC_BARRIER = 3, + VAIDYA_BARRIER = 4 +}; + +template +struct AssertBarrierFalseType : std::false_type {}; + template struct AssertFalseType : std::false_type {}; +template +auto get_max_step(VT const& Ad, VT const& b_Ax) +{ + using NT = typename VT::Scalar; + const int m = Ad.size(); + NT max_element = std::numeric_limits::lowest(), max_element_temp; + for (int i = 0; i < m; i++) { + max_element_temp = Ad.coeff(i) / b_Ax.coeff(i); + if (max_element_temp > max_element) { + max_element = max_element_temp; + } + } + + return NT(1) / max_element; +} + template inline static auto initialize_chol(MT const& mat) @@ -29,7 +57,7 @@ initialize_chol(MT const& mat) if constexpr (std::is_same::value) { return std::make_unique>(); - } else if constexpr (std::is_same::value) + } else if constexpr (std::is_base_of, MT >::value) { auto llt = std::make_unique>(); llt->analyzePattern(mat); @@ -50,7 +78,7 @@ initialize_chol(MT const& A_trans, MT const& A) if constexpr (std::is_same::value) { return std::make_unique>(); - } else if constexpr (std::is_same::value) + } else if constexpr (std::is_base_of, MT >::value) { MT mat = A_trans * A; return initialize_chol(mat); @@ -71,7 +99,7 @@ inline static VT solve_vec(std::unique_ptr const& llt, { llt->compute(H); return llt->solve(b); - } else if constexpr (std::is_same::value) + } else if constexpr (std::is_base_of, MT >::value) { llt->factorize(H); return llt->solve(b); @@ -94,7 +122,7 @@ solve_mat(std::unique_ptr const& llt, llt->compute(H); logdetE = llt->matrixL().toDenseMatrix().diagonal().array().log().sum(); return llt->solve(mat); - } else if constexpr (std::is_same::value) + } else if constexpr (std::is_base_of, MT >::value) { llt->factorize(H); logdetE = llt->matrixL().nestedExpression().diagonal().array().log().sum(); @@ -115,7 +143,7 @@ inline static void update_Atrans_Diag_A(MT &H, MT const& A_trans, if constexpr (std::is_same::value) { H.noalias() = A_trans * D * A; - } else if constexpr (std::is_same::value) + } else if constexpr (std::is_base_of, MT >::value) { H = A_trans * D * A; } else @@ -133,7 +161,7 @@ inline static void update_Diag_A(MT &H, diag_MT const& D, MT const& A) if constexpr (std::is_same::value) { H.noalias() = D * A; - } else if constexpr (std::is_same::value) + } else if constexpr (std::is_base_of, MT >::value) { H = D * A; } else @@ -151,7 +179,7 @@ inline static void update_A_Diag(MT &H, MT const& A, diag_MT const& D) if constexpr (std::is_same::value) { H.noalias() = A * D; - } else if constexpr (std::is_same::value) + } else if constexpr (std::is_base_of, MT >::value) { H = A * D; } else @@ -170,7 +198,7 @@ get_mat_prod_op(MT const& E) if constexpr (std::is_same::value) { return std::make_unique>(E); - } else if constexpr (std::is_same::value) + } else if constexpr (std::is_base_of, MT >::value) { return std::make_unique>(E); } else @@ -221,7 +249,7 @@ init_Bmat(MT &B, int const n, MT const& A_trans, MT const& A) if constexpr (std::is_same::value) { B.resize(n+1, n+1); - } else if constexpr (std::is_same::value) + } else if constexpr (std::is_base_of, MT >::value) { // Initialize the structure of matrix B typedef Eigen::Triplet triplet; @@ -271,7 +299,7 @@ update_Bmat(MT &B, VT const& AtDe, VT const& d, B.block(n, 0, 1, n).noalias() = AtDe.transpose(); B(n, n) = d.sum(); B.noalias() += 1e-14 * MT::Identity(n + 1, n + 1); - } else if constexpr (std::is_same::value) + } else if constexpr (std::is_base_of, MT >::value) { MT AtD_A = AtD * A; int k = 0; @@ -312,5 +340,91 @@ update_Bmat(MT &B, VT const& AtDe, VT const& d, } } +template +std::tuple init_step() +{ + if constexpr (BarrierType == EllipsoidType::LOG_BARRIER) + { + return {NT(1), NT(0.99)}; + } else if constexpr (BarrierType == EllipsoidType::VOLUMETRIC_BARRIER || + BarrierType == EllipsoidType::VAIDYA_BARRIER) + { + return {NT(0.5), NT(0.4)}; + } else { + static_assert(AssertBarrierFalseType::value, + "Barrier type is not supported."); + } +} + +template +void get_barrier_hessian_grad(MT const& A, MT const& A_trans, VT const& b, + VT const& x, VT const& Ax, llt_type const& llt, + MT &H, VT &grad, VT &b_Ax, NT &obj_val) +{ + b_Ax.noalias() = b - Ax; + VT s = b_Ax.cwiseInverse(); + VT s_sq = s.cwiseProduct(s); + VT sigma; + // Hessian of the log-barrier function + update_Atrans_Diag_A(H, A_trans, A, s_sq.asDiagonal()); + + if constexpr (BarrierType == EllipsoidType::VOLUMETRIC_BARRIER || + BarrierType == EllipsoidType::VAIDYA_BARRIER) + { + // Computing sigma(x)_i = (a_i^T H^{-1} a_i) / (b_i - a_i^Tx)^2 + MT_dense HA = solve_mat(llt, H, A_trans, obj_val); + MT_dense aiHai = HA.transpose().cwiseProduct(A); + sigma = (aiHai.rowwise().sum()).cwiseProduct(s_sq); + } + + if constexpr (BarrierType == EllipsoidType::LOG_BARRIER) + { + grad.noalias() = A_trans * s; + } else if constexpr (BarrierType == EllipsoidType::VOLUMETRIC_BARRIER) + { + // Gradient of the volumetric barrier function + grad.noalias() = A_trans * (s.cwiseProduct(sigma)); + // Hessian of the volumetric barrier function + update_Atrans_Diag_A(H, A_trans, A, s_sq.cwiseProduct(sigma).asDiagonal()); + } else if constexpr (BarrierType == EllipsoidType::VAIDYA_BARRIER) + { + const int m = b.size(), d = x.size(); + NT const d_m = NT(d) / NT(m); + // Weighted gradient of the log barrier function + grad.noalias() = A_trans * s; + grad *= d_m; + // Add the gradient of the volumetric function + grad.noalias() += A_trans * (s.cwiseProduct(sigma)); + // Weighted Hessian of the log barrier function + H *= d_m; + // Add the Hessian of the volumetric function + MT Hvol(d, d); + update_Atrans_Diag_A(Hvol, A_trans, A, s_sq.cwiseProduct(sigma).asDiagonal()); + H += Hvol; + obj_val -= s.array().log().sum(); + } else { + static_assert(AssertBarrierFalseType::value, + "Barrier type is not supported."); + } +} + +template +void get_step_next_iteration(NT const obj_val_prev, NT const obj_val, + NT const tol_obj, NT &step_iter) +{ + if constexpr (BarrierType == EllipsoidType::LOG_BARRIER) + { + step_iter = NT(1); + } else if constexpr (BarrierType == EllipsoidType::VOLUMETRIC_BARRIER) + { + step_iter *= (obj_val_prev <= obj_val - tol_obj) ? NT(0.9) : NT(0.999); + } else if constexpr (BarrierType == EllipsoidType::VAIDYA_BARRIER) + { + step_iter *= NT(0.999); + } else { + static_assert(AssertBarrierFalseType::value, + "Barrier type is not supported."); + } +} -#endif // MAT_COMPUTATIONAL_OPERATORS_H +#endif // ROUNDING_UTIL_FUNCTIONS_HPP \ No newline at end of file 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) +NT get_next_gaussian(Polytope& P, + Point &p, + NT const& a, + const unsigned int &N, + const NT &ratio, + const NT &C, + const unsigned int& walk_length, + RandomNumberGenerator& rng, + Grad& g, + Func& f, + crhmc_walk_params& parameters, + CrhmcProblem& problem, + CRHMCWalkType& crhmc_walk) +{ + + NT last_a = a; + NT last_ratio = 0.1; + //k is needed for the computation of the next variance a_{i+1} = a_i * (1-1/d)^k + NT k = 1.0; + const NT tol = 0.00001; + bool done=false; + std::vector fn(N,NT(0.0)); + std::list randPoints; + typedef typename std::vector::iterator viterator; + + //sample N points + PushBackWalkPolicy push_back_policy; + bool raw_output = false; + + typedef CrhmcRandomPointGenerator CRHMCRandomPointGenerator; + + CRHMCRandomPointGenerator::apply(problem, p, N, walk_length, randPoints, + push_back_policy, rng, g, f, parameters, crhmc_walk, simdLen, raw_output); + + while (!done) + { + NT new_a = last_a * std::pow(ratio,k); + + auto fnit = fn.begin(); + for (auto pit=randPoints.begin(); pit!=randPoints.end(); ++pit, fnit++) + { + *fnit = eval_exp(*pit, new_a)/eval_exp(*pit, last_a); + } + std::pair mv = get_mean_variance(fn); + + // Compute a_{i+1} + if (mv.second/(mv.first * mv.first)>=C || mv.first/last_ratio<1.0+tol) + { + if (k != 1.0) + { + k = k / 2; + } + done = true; + } + else { + k = 2 * k; + } + last_ratio = mv.first; + } + return last_a * std::pow(ratio, k); +} + +// Compute the sequence of spherical gaussians +template +< + int simdLen, + typename Polytope, + typename NT, + typename RandomNumberGenerator +> +void compute_annealing_schedule(Polytope& P, + NT const& ratio, + NT const& C, + NT const& frac, + unsigned int const& N, + unsigned int const& walk_length, + NT const& chebychev_radius, + NT const& error, + std::vector& a_vals, + RandomNumberGenerator& rng) +{ + typedef typename Polytope::PointType Point; + typedef typename Polytope::MT MT; + + typedef typename GaussianFunctor::FunctionFunctor Func; + typedef typename GaussianFunctor::GradientFunctor Grad; + typedef typename GaussianFunctor::HessianFunctor Hess; + typedef typename GaussianFunctor::parameters func_params; + + typedef crhmc_input Input; + typedef crhmc_problem CrhmcProblem; + + typedef ImplicitMidpointODESolver Solver; + + typedef typename CRHMCWalk::template Walk + < + Point, + CrhmcProblem, + RandomNumberGenerator, + Grad, + Func, + Solver + > CRHMCWalkType; + + typedef typename CRHMCWalk::template parameters + < + NT, + Grad + > crhmc_walk_params; + + typedef CrhmcRandomPointGenerator CRHMCRandomPointGenerator; + + + // Compute the first gaussian + // This uses the function from the standard volume_cooling_gaussians.hpp + get_first_gaussian(P, frac, chebychev_radius, error, a_vals); + NT a_stop = 0.0; + const NT tol = 0.001; + unsigned int it = 0; + unsigned int n = P.dimension(); + const unsigned int totalSteps = ((int)150/((1.0 - frac) * error))+1; + + if (a_vals[0](P, f, g, h); + + typedef crhmc_problem CrhmcProblem; + CrhmcProblem problem = CrhmcProblem(input); + + Point p = Point(problem.center); + + if(problem.terminate){return;} + + problem.options.simdLen = simdLen; + crhmc_walk_params params(input.df, p.dimension(), problem.options); + + if (input.df.params.eta > 0) { + params.eta = input.df.params.eta; + } + + int dim; + dim = p.dimension(); + + //create the walk object for this problem + CRHMCWalkType walk = CRHMCWalkType(problem, p, input.df, input.f, params); + + // Compute the next gaussian + NT next_a = get_next_gaussian + (P, p, a_vals[it], N, ratio, C, walk_length, rng, g, f, params, problem, walk); + +#ifdef VOLESTI_DEBUG + std::cout<<"Next Gaussian " << next_a <0 && curr_fn/curr_its>(1.0+tol)) + { + a_vals.push_back(next_a); + it++; + } else if (next_a <= 0) + { + a_vals.push_back(a_stop); + it++; + break; + } else { + a_vals[it] = a_stop; + break; + } + } +#ifdef VOLESTI_DEBUG + std::cout<<"first gaussian after WHILE "<< a_vals[0] < +double volume_cooling_gaussians(Polytope& Pin, + RandomNumberGenerator& rng, + double const& error = 0.1, + unsigned int const& walk_length = 1) +{ + typedef typename Polytope::PointType Point; + typedef typename Point::FT NT; + typedef typename Polytope::VT VT; + typedef typename Polytope::MT MT; + typedef typename GaussianFunctor::FunctionFunctor Func; + typedef typename GaussianFunctor::GradientFunctor Grad; + typedef typename GaussianFunctor::HessianFunctor Hess; + typedef typename GaussianFunctor::parameters func_params; + + typedef crhmc_input Input; + typedef crhmc_problem CrhmcProblem; + + typedef ImplicitMidpointODESolver Solver; + + typedef typename CRHMCWalk::template Walk + < + Point, + CrhmcProblem, + RandomNumberGenerator, + Grad, + Func, + Solver + > CRHMCWalkType; + + typedef typename CRHMCWalk::template parameters + < + NT, + Grad + > crhmc_walk_params; + + typedef CrhmcRandomPointGenerator CRHMCRandomPointGenerator; + + auto P(Pin); //copy and work with P because we are going to shift + unsigned int n = P.dimension(); + unsigned int m = P.num_of_hyperplanes(); + gaussian_annealing_parameters parameters(P.dimension()); + + // Consider Chebychev center as an internal point + auto InnerBall = P.ComputeInnerBall(); + if (InnerBall.second < 0.0) return -1.0; + + Point c = InnerBall.first; + NT radius = InnerBall.second; + + // Move the chebychev center to the origin and apply the same shifting to the polytope + P.shift(c.getCoefficients()); + + // Computing the sequence of gaussians +#ifdef VOLESTI_DEBUG + std::cout<<"\n\nComputing annealing...\n"< a_vals; + NT ratio = parameters.ratio; + NT C = parameters.C; + unsigned int N = parameters.N; + + compute_annealing_schedule(P, ratio, C, parameters.frac, N, walk_length, radius, error, a_vals, rng); + +#ifdef VOLESTI_DEBUG + std::cout<<"All the variances of schedule_annealing computed in = " + << (double)clock()/(double)CLOCKS_PER_SEC-tstart2<<" sec"< last_W2(W,0); + std::vector fn(mm,0); + std::vector its(mm,0); + VT lamdas; + lamdas.setZero(m); + NT vol = std::pow(M_PI/a_vals[0], (NT(n))/2.0); + unsigned int i=0; + + typedef typename std::vector::iterator viterator; + viterator itsIt = its.begin(); + viterator avalsIt = a_vals.begin(); + viterator minmaxIt; + + +#ifdef VOLESTI_DEBUG + std::cout<<"volume of the first gaussian = "<::min(); + NT max_val = std::numeric_limits::max(); + unsigned int min_index = W-1; + unsigned int max_index = W-1; + unsigned int index = 0; + unsigned int min_steps = 0; + std::vector last_W = last_W2; + + // Set the radius for the ball walk + //creating the walk object + int dimension = P.dimension(); + func_params f_params = func_params(Point(dimension), *avalsIt, 1); + + Func f(f_params); + Grad g(f_params); + Hess h(f_params); + + //create the crhmc problem + Input input = convert2crhmc_input(P, f, g, h); + + CrhmcProblem problem = CrhmcProblem(input); + + Point p = Point(problem.center); + + if(problem.terminate){return 0;} + + problem.options.simdLen=simdLen; + crhmc_walk_params params(input.df, p.dimension(), problem.options); + + if (input.df.params.eta > 0) { + params.eta = input.df.params.eta; + } + + CRHMCWalkType walk = CRHMCWalkType(problem, p, input.df, input.f, params); + + while (!done || (*itsIt)= max_val) + { + max_val = val; + max_index = index; + } else if (max_index == index) + { + minmaxIt = std::max_element(last_W.begin(), last_W.end()); + max_val = *minmaxIt; + max_index = std::distance(last_W.begin(), minmaxIt); + } + + if ( (max_val-min_val)/max_val <= curr_eps/2.0 ) + { + done=true; + } + + index = index%W + 1; + if (index == W) index = 0; + } +#ifdef VOLESTI_DEBUG + std::cout << "ratio " << i << " = " << (*fnIt) / (*itsIt) + << " N_" << i << " = " << *itsIt << std::endl; +#endif + vol *= ((*fnIt) / (*itsIt)); + } + +#ifdef VOLESTI_DEBUG + NT sum_of_steps = 0.0; + for(viterator it = its.begin(); it != its.end(); ++it) { + sum_of_steps += *it; + } + auto steps= int(sum_of_steps); + std::cout<<"\nTotal number of steps = "< +void burnIn(Point &p, + const unsigned int& walk_length, + RandomNumberGenerator& rng, + Grad& g, + Func& f, + crhmc_walk_params& parameters, + CrhmcProblem& problem, + CRHMCWalkType& crhmc_walk, + NT burnIn_sample) +{ + std::list randPoints; + + // burnIn + PushBackWalkPolicy push_back_policy; + bool raw_output = false; + typedef CrhmcRandomPointGenerator CRHMCRandomPointGenerator; + + CRHMCRandomPointGenerator::apply(problem, p, burnIn_sample, walk_length, randPoints, + push_back_policy, rng, g, f, parameters, crhmc_walk, simdLen, raw_output); + +} + + +template +< + int simdLen, + typename Polytope, + typename Point, + typename NT, + typename MT, + typename RandomNumberGenerator +> +NT get_next_gaussian(Polytope& P, + Point& p, + NT const& a, + const unsigned int &N, + const NT &ratio, + const NT &C, + const unsigned int& walk_length, + RandomNumberGenerator& rng, + MT const& inv_covariance_matrix) +{ + typedef typename NonSphericalGaussianFunctor::FunctionFunctor Func; + typedef typename NonSphericalGaussianFunctor::GradientFunctor Grad; + typedef typename NonSphericalGaussianFunctor::HessianFunctor Hess; + typedef typename NonSphericalGaussianFunctor::parameters func_params; + + typedef crhmc_input> Input; + + typedef crhmc_problem CrhmcProblem; + + typedef ImplicitMidpointODESolver Solver; + + typedef typename CRHMCWalk::template Walk + < + Point, + CrhmcProblem, + RandomNumberGenerator, + Grad, + Func, + Solver + > CRHMCWalkType; + + typedef typename CRHMCWalk::template parameters + < + NT, + Grad + > crhmc_walk_params; + + typedef CrhmcRandomPointGenerator CRHMCRandomPointGenerator; + + // Create the CRHMC problem for this variance + int dimension = P.dimension(); + func_params f_params = func_params(Point(dimension), a, -1, inv_covariance_matrix); + + Func f(f_params); + Grad g(f_params); + Hess h(f_params); + ZeroFunctor zerof; + + Input input = convert2crhmc_input>(P, f, g, zerof); + + typedef crhmc_problem CrhmcProblem; + + CrhmcProblem problem = CrhmcProblem(input); + + if(problem.terminate) { return 0;} + + problem.options.simdLen = simdLen; + + crhmc_walk_params params(input.df, p.dimension(), problem.options); + + // Create the walk object for this problem + CRHMCWalkType walk = CRHMCWalkType(problem, p, input.df, input.f, params); + + burnIn( + p, walk_length, rng, g, f, params, problem, walk, 2.5 * N); + + NT last_a = a; + NT last_ratio = 0.1; + NT k = 1.0; + const NT tol = 0.00001; + bool done = false; + std::vector fn(N, NT(0.0)); + std::list randPoints; + + // sample N points + PushBackWalkPolicy push_back_policy; + bool raw_output = false; + typedef CrhmcRandomPointGenerator CRHMCRandomPointGenerator; + + CRHMCRandomPointGenerator::apply(problem, p, N, walk_length, randPoints, + push_back_policy, rng, g, f, params, walk, simdLen, raw_output); + + while (!done) { + NT new_a = last_a * std::pow(ratio, k); + + auto fnit = fn.begin(); + for (auto pit = randPoints.begin(); pit != randPoints.end(); ++pit, fnit++) { + *fnit = eval_exp(*pit, inv_covariance_matrix, new_a, last_a); + } + + std::pair mv = get_mean_variance(fn); + + // Compute a_{i+1} + if (mv.second / (mv.first * mv.first) >= C || mv.first / last_ratio < 1.0 + tol) { + if (k != 1.0) { + k = k / 2; + } + done = true; + } else { + k = 2 * k; + } + + last_ratio = mv.first; + } + + // Return the new a value as a scalar + return last_a * std::pow(ratio, k); +} + +// Compute the sequence of non spherical gaussians +template< + int simdLen, + typename Point, + typename Polytope, + typename NT, + typename MT, + typename RandomNumberGenerator +> +void compute_annealing_schedule(Polytope Pin_copy, + Polytope P_copy, + NT const& ratio, + NT const& C, + NT const& frac, + unsigned int const& N, + unsigned int const& walk_length, + NT const& error, + std::vector& a_vals, + MT const& inv_covariance_matrix, + RandomNumberGenerator& rng) +{ + + typedef typename NonSphericalGaussianFunctor::FunctionFunctor Func; + typedef typename NonSphericalGaussianFunctor::GradientFunctor Grad; + typedef typename NonSphericalGaussianFunctor::HessianFunctor Hess; + typedef typename NonSphericalGaussianFunctor::parameters func_params; + + typedef crhmc_input> Input; + + typedef crhmc_problem CrhmcProblem; + + typedef ImplicitMidpointODESolver Solver; + + typedef typename CRHMCWalk::template Walk + < + Point, + CrhmcProblem, + RandomNumberGenerator, + Grad, + Func, + Solver + > CRHMCWalkType; + + typedef typename CRHMCWalk::template parameters + < + NT, + Grad + > crhmc_walk_params; + + typedef CrhmcRandomPointGenerator CRHMCRandomPointGenerator; + + // compute the first variance + // for this we need P + auto ball = P_copy.ComputeInnerBall(); + P_copy.shift(ball.first.getCoefficients()); // when using max_ellipsoid for rounding this center is the origin, but if we use other covariances this is different than the origin + get_first_gaussian(P_copy, frac, ball.second, error, a_vals); // the function included from volume_cooling_gaussians.hpp (spherical gaussians) +#ifdef VOLESTI_DEBUG + std::cout << "first gaussian computed " << a_vals[0] << std::endl; +#endif + + //for the rest we need Pin + NT a_stop = 0.0; + const NT tol = 0.001; + unsigned int it = 0; + unsigned int n = Pin_copy.dimension(); + const unsigned int totalSteps = ((int)150/((1.0 - frac) * error)) + 1; + + if (a_vals[0] initial_zerof; + + Input initial_input = convert2crhmc_input>(Pin_copy, initial_f, initial_g, initial_zerof); + CrhmcProblem initial_problem = CrhmcProblem(initial_input); + + Point initial_p = Point(initial_problem.center); + initial_problem.options.simdLen = simdLen; + crhmc_walk_params initial_params(initial_input.df, initial_p.dimension(), initial_problem.options); + CRHMCWalkType initial_walk = CRHMCWalkType(initial_problem, initial_p, initial_input.df, initial_input.f, initial_params); + + burnIn( + initial_p, walk_length, rng, initial_g, initial_f, initial_params, initial_problem, initial_walk, 2.5 * N); + + //fix eta + initial_eta = initial_walk.get_current_eta(); + + //fix point + start_point = initial_problem.center; + + while (true) { + + // Compute the next gaussian + NT next_a = get_next_gaussian( + Pin_copy, start_point, a_vals[it], N, ratio, C, walk_length, rng, inv_covariance_matrix); + + // Decide if this is the last one + NT curr_fn = 0; + NT curr_its = 0; + auto steps = totalSteps; + + //TODO: potential problem creation and preprocessing optimization + + // Create the CRHMC problem for this variance + int dimension = Pin_copy.dimension(); + func_params f_params = func_params(Point(dimension), a_vals[it], -1, inv_covariance_matrix); + + Func f(f_params); + Grad g(f_params); + ZeroFunctor zerof; + + Input input = convert2crhmc_input>(Pin_copy, f, g, zerof); + + typedef crhmc_problem CrhmcProblem; + + CrhmcProblem problem = CrhmcProblem(input); + + if(problem.terminate) { return; } + + problem.options.simdLen = simdLen; + + crhmc_walk_params params(input.df, start_point.dimension(), problem.options); + + //eta initialization with the previous walk eta + problem.options.etaInitialize = false; + params.eta = initial_eta; + + // Create the walk object for this problem + CRHMCWalkType walk = CRHMCWalkType(problem, start_point, input.df, input.f, params); + + // Compute some ratios to decide if this is the last gaussian + for (unsigned int j = 0; j < steps; j++) + { + walk.apply(rng, walk_length); + start_point = walk.getPoint(); + curr_its += 1.0; + curr_fn += eval_exp(start_point, inv_covariance_matrix, next_a, a_vals[it]); + } + + //restore the new eta and start point, by looking at the walk after its operations + initial_eta = walk.get_current_eta(); + + // Remove the last gaussian. + // Set the last a_i equal to zero + if (next_a > 0 && curr_fn / curr_its > (1.0 + tol)) + { + a_vals.push_back(next_a); + it++; + } else if (next_a <= 0) + { + a_vals.push_back(a_stop); + it++; + break; + } else { + a_vals[it] = a_stop; + break; + } + } +} + + +template +struct non_gaussian_annealing_parameters +{ + non_gaussian_annealing_parameters(unsigned int d) + : frac(0.1) + , ratio(NT(1) - NT(1) / NT(d)) + , C(NT(2)) + , N(500 * ((int) C) + ((int) (d * d))) + , W(6 * d * d + 800) + {} + + NT frac; + NT ratio; + NT C; + unsigned int N; + unsigned int W; +}; + + +template +< + typename Polytope, + typename RandomNumberGenerator, + typename WalkTypePolicy = CRHMCWalk, + int simdLen = 8 +> +double non_spherical_crhmc_volume_cooling_gaussians(Polytope& Pin, + RandomNumberGenerator& rng, + double const& error = 0.1, + unsigned int const& walk_length = 1) +{ + typedef typename Polytope::PointType Point; + typedef typename Point::FT NT; + typedef typename Polytope::VT VT; + typedef typename Polytope::MT MT; + typedef typename NonSphericalGaussianFunctor::FunctionFunctor Func; + typedef typename NonSphericalGaussianFunctor::GradientFunctor Grad; + typedef typename NonSphericalGaussianFunctor::parameters func_params; + + typedef crhmc_input> Input; + typedef crhmc_problem CrhmcProblem; + + typedef ImplicitMidpointODESolver Solver; + + typedef typename WalkTypePolicy::template Walk + < + Point, + CrhmcProblem, + RandomNumberGenerator, + Grad, + Func, + Solver + > CRHMCWalkType; + + typedef typename WalkTypePolicy::template parameters + < + NT, + Grad + > crhmc_walk_params; + + typedef HPolytope HPOLYTOPE; + + HPOLYTOPE P(Pin.dimension(), Pin.get_mat(), Pin.get_vec()); + HPOLYTOPE newPin(Pin.dimension(), Pin.get_mat(), Pin.get_vec()); + unsigned int n = P.dimension(); + unsigned int m = P.num_of_hyperplanes(); + + //compute inscribed ellipsoid + NT tol = std::pow(10, -6.0), reg = std::pow(10, -4.0); + unsigned int maxiter = 100; + P.normalize(); + VT x0 = compute_feasible_point(P.get_mat(), P.get_vec()); + auto ellipsoid_result = compute_inscribed_ellipsoid(P.get_mat(), P.get_vec(), x0, maxiter, tol, reg); + + // extract the covariance matrix and the center of the ellipsoid + MT inv_covariance_matrix = std::get<0>(ellipsoid_result); //this is the covariance to use in the telescopic product + VT center = std::get<1>(ellipsoid_result); + MT covariance_matrix = inv_covariance_matrix.inverse(); + + newPin.shift(center); //we shift the initial polytope so that the origin is the center of the gaussians + P.shift(center); + + // we apply the rounding transformation + Eigen::LLT lltOfA(covariance_matrix); + auto L = lltOfA.matrixL(); + P.linear_transformIt(L); + + + // Initialize the gaussian_annealing_parameters struct + non_gaussian_annealing_parameters parameters(P.dimension()); + + // Computing the sequence of gaussians +#ifdef VOLESTI_DEBUG + std::cout<<"\n\nComputing annealing...\n"< a_vals; + NT ratio = parameters.ratio; + NT C = parameters.C; + unsigned int N = parameters.N; + + compute_annealing_schedule(newPin, P, ratio, C, parameters.frac, N, walk_length, error, a_vals, inv_covariance_matrix, rng); + +#ifdef VOLESTI_DEBUG + std::cout<<"All the variances of schedule_annealing computed in = " + << (double)clock()/(double)CLOCKS_PER_SEC-tstart2<<" sec"< last_W2(W,0); + std::vector fn(mm,0); + std::vector its(mm,0); + VT lamdas; + lamdas.setZero(m); + + MT scaled_matrix = a_vals[0] * inv_covariance_matrix; + NT det_scaled = scaled_matrix.determinant(); + NT vol = std::pow(M_PI, n / 2.0) / std::sqrt(det_scaled); + + unsigned int i=0; + + typedef typename std::vector::iterator viterator; + viterator itsIt = its.begin(); + auto avalsIt = a_vals.begin(); + viterator minmaxIt; + + +#ifdef VOLESTI_DEBUG + std::cout<<"volume of the first gaussian = "<::min(); + NT max_val = std::numeric_limits::max(); + unsigned int min_index = W-1; + unsigned int max_index = W-1; + unsigned int index = 0; + unsigned int min_steps = 0; + std::vector last_W = last_W2; + + // Set the radius for the ball walk + //creating the walk object + int dimension = newPin.dimension(); + func_params f_params = func_params(Point(dimension), *avalsIt, -1, inv_covariance_matrix); + + Func f(f_params); + Grad g(f_params); + + ZeroFunctor zerof; + + //create the crhmc problem + Input input = convert2crhmc_input>(newPin, f, g, zerof); + + typedef crhmc_problem CrhmcProblem; + CrhmcProblem problem = CrhmcProblem(input); + Point p = problem.center; + + if(problem.terminate){return 0;} + problem.options.simdLen=simdLen; + + //create the walk and do the burnIn + crhmc_walk_params params(input.df, p.dimension(), problem.options); + CRHMCWalkType walk = CRHMCWalkType(problem, p, input.df, input.f, params); + + burnIn( + p, walk_length, rng, g, f, params, problem, walk, 2.5 * N); + + while (!done || (*itsIt) < min_steps) + { + walk.apply(rng, walk_length); + p = walk.getPoint(); + *itsIt = *itsIt + 1.0; + + *fnIt += eval_exp(p, inv_covariance_matrix, *(avalsIt+1), *avalsIt); + + NT val = (*fnIt) / (*itsIt); + + last_W[index] = val; + if (val <= min_val) + { + min_val = val; + min_index = index; + } else if (min_index == index) + { + minmaxIt = std::min_element(last_W.begin(), last_W.end()); + min_val = *minmaxIt; + min_index = std::distance(last_W.begin(), minmaxIt); + } + + if (val >= max_val) + { + max_val = val; + max_index = index; + } else if (max_index == index) + { + minmaxIt = std::max_element(last_W.begin(), last_W.end()); + max_val = *minmaxIt; + max_index = std::distance(last_W.begin(), minmaxIt); + } + + if ((max_val - min_val) / max_val <= curr_eps / 2.0) + { + done = true; + } + + index = index % W + 1; + if (index == W) index = 0; + } +#ifdef VOLESTI_DEBUG + std::cout << "ratio " << i << " = " << (*fnIt) / (*itsIt) + << " N_" << i << " = " << *itsIt << std::endl; +#endif + vol *= ((*fnIt) / (*itsIt)); + } + +#ifdef VOLESTI_DEBUG + NT sum_of_steps = 0.0; + for(viterator it = its.begin(); it != its.end(); ++it) { + sum_of_steps += *it; + } + auto steps= int(sum_of_steps); + std::cout<<"\nTotal number of steps = "<