From 9d2eb4c60c3df5ff985594fc4c62eb72df3ad266 Mon Sep 17 00:00:00 2001 From: lucaperju <104904648+lucaperju@users.noreply.github.com> Date: Mon, 22 Jul 2024 17:42:18 +0300 Subject: [PATCH] Gaussian Accelerated Billiard Walk (#323) * New gaussian accelerated billiard walk * minor changes for PR * add space * Fix indentation, simplify equation * minor optimizations * added normalized flag, made indentation consistent, changed file name * remove line * add copyright comments * added is_normalized function to all types of bodies * remove old file * fix outdated example --- .../sampler.cpp | 16 +- include/convex_bodies/ballintersectconvex.h | 6 + .../correlation_spectrahedron.hpp | 6 + include/convex_bodies/hpolytope.h | 25 ++- include/convex_bodies/orderpolytope.h | 4 + .../spectrahedra/spectrahedron.h | 7 + include/convex_bodies/vpolyintersectvpoly.h | 4 + include/convex_bodies/vpolytope.h | 4 + include/convex_bodies/zonoIntersecthpoly.h | 6 + include/convex_bodies/zpolytope.h | 5 + .../gaussian_accelerated_billiard_walk.hpp | 187 ++++++++++-------- .../uniform_accelerated_billiard_walk.hpp | 38 ++-- ...orm_accelerated_billiard_walk_parallel.hpp | 23 ++- include/sampling/random_point_generators.hpp | 4 +- test/CMakeLists.txt | 1 + test/sampling_test.cpp | 53 +++++ 16 files changed, 274 insertions(+), 115 deletions(-) diff --git a/examples/sampling-hpolytope-with-billiard-walks/sampler.cpp b/examples/sampling-hpolytope-with-billiard-walks/sampler.cpp index 803d38dc3..efdc4a664 100644 --- a/examples/sampling-hpolytope-with-billiard-walks/sampler.cpp +++ b/examples/sampling-hpolytope-with-billiard-walks/sampler.cpp @@ -7,6 +7,7 @@ #include "sampling/random_point_generators.hpp" #include "random_walks/random_walks.hpp" #include "preprocess/max_inscribed_ellipsoid.hpp" +#include "preprocess/inscribed_ellipsoid_rounding.hpp" #include "convex_bodies/ellipsoid.h" #include "convex_bodies/hpolytope.h" @@ -69,23 +70,18 @@ void sample_using_gaussian_billiard_walk(HPOLYTOPE& HP, RNGType& rng, unsigned i Point q(HP.dimension()); // origin // ----------- Get inscribed ellipsoid -------------------------------- - typedef Ellipsoid EllipsoidType; unsigned int max_iter = 150; NT tol = std::pow(10, -6.0), reg = std::pow(10, -4.0); VT x0 = q.getCoefficients(); - MT E; VT center; bool converged; - std::tie(E, center, converged) = max_inscribed_ellipsoid(HP.get_mat(), - HP.get_vec(), x0, max_iter, tol, reg); - - if (!converged) // not converged - throw std::runtime_error("max_inscribed_ellipsoid not converged"); - - EllipsoidType inscribed_ellipsoid(E); + std::tuple ellipsoid = compute_inscribed_ellipsoid + (HP.get_mat(), HP.get_vec(), x0, max_iter, tol, reg); + + const MT E = get<0>(ellipsoid); // -------------------------------------------------------------------- - Generator::apply(HP, q, inscribed_ellipsoid, num_points, walk_len, + Generator::apply(HP, q, E, num_points, walk_len, randPoints, push_back_policy, rng); write_to_file(filename, randPoints); } diff --git a/include/convex_bodies/ballintersectconvex.h b/include/convex_bodies/ballintersectconvex.h index b61ae5324..992aae1a1 100644 --- a/include/convex_bodies/ballintersectconvex.h +++ b/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/include/convex_bodies/correlation_matrices/correlation_spectrahedron.hpp b/include/convex_bodies/correlation_matrices/correlation_spectrahedron.hpp index a615b33c0..92c38dac6 100755 --- a/include/convex_bodies/correlation_matrices/correlation_spectrahedron.hpp +++ b/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/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h index d0e894faa..b01b6307c 100644 --- a/include/convex_bodies/hpolytope.h +++ b/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 @@ -56,6 +57,7 @@ class HPolytope { 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 public: //TODO: the default implementation of the Big3 should be ok. Recheck. @@ -68,7 +70,7 @@ class HPolytope { // Copy constructor HPolytope(HPolytope const& p) : - _d{p._d}, A{p.A}, b{p.b}, _inner_ball{p._inner_ball} + _d{p._d}, A{p.A}, b{p.b}, _inner_ball{p._inner_ball}, normalized{p.normalized} { } @@ -172,11 +174,17 @@ class HPolytope { return b; } + bool is_normalized () + { + return normalized; + } + // change the matrix A void set_mat(MT const& A2) { A = A2; + normalized = false; } @@ -541,7 +549,6 @@ class HPolytope { NT lamda = 0; VT sum_nom; - unsigned int j; int m = num_of_hyperplanes(), facet; Ar.noalias() += lambda_prev*Av; @@ -822,6 +829,7 @@ class HPolytope { void linear_transformIt(MT const& T) { A = A * T; + normalized = false; } @@ -866,6 +874,7 @@ class HPolytope { A.row(i) = A.row(i) / row_norm; b(i) = b(i) / row_norm; } + normalized = true; } void compute_reflection(Point& v, Point const&, int const& facet) const @@ -909,6 +918,18 @@ class HPolytope { v += a; } + template + NT compute_reflection(Point &v, const Point &, MT 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&){} template diff --git a/include/convex_bodies/orderpolytope.h b/include/convex_bodies/orderpolytope.h index 79f41e161..12489f247 100644 --- a/include/convex_bodies/orderpolytope.h +++ b/include/convex_bodies/orderpolytope.h @@ -100,6 +100,10 @@ class OrderPolytope { return b; } + bool is_normalized () + { + return _normalized; + } // print polytope in Ax <= b format void print() const diff --git a/include/convex_bodies/spectrahedra/spectrahedron.h b/include/convex_bodies/spectrahedra/spectrahedron.h index e8f341382..37b997ea1 100644 --- a/include/convex_bodies/spectrahedra/spectrahedron.h +++ b/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/include/convex_bodies/vpolyintersectvpoly.h b/include/convex_bodies/vpolyintersectvpoly.h index a6dd97caf..ec2a7c5dd 100644 --- a/include/convex_bodies/vpolyintersectvpoly.h +++ b/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/include/convex_bodies/vpolytope.h b/include/convex_bodies/vpolytope.h index 82ca66891..9100756ef 100644 --- a/include/convex_bodies/vpolytope.h +++ b/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/include/convex_bodies/zonoIntersecthpoly.h b/include/convex_bodies/zonoIntersecthpoly.h index bbad43a68..e356db8db 100644 --- a/include/convex_bodies/zonoIntersecthpoly.h +++ b/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/include/convex_bodies/zpolytope.h b/include/convex_bodies/zpolytope.h index 1ad8868ea..68ce0cf0e 100644 --- a/include/convex_bodies/zpolytope.h +++ b/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/include/random_walks/gaussian_accelerated_billiard_walk.hpp b/include/random_walks/gaussian_accelerated_billiard_walk.hpp index fbf43f443..ea059305a 100644 --- a/include/random_walks/gaussian_accelerated_billiard_walk.hpp +++ b/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,67 @@ struct GaussianAcceleratedBilliardWalk template - < - typename Polytope, - typename RandomNumberGenerator - > + < + typename Polytope, + typename RandomNumberGenerator + > struct Walk { typedef typename Polytope::PointType Point; - // typedef typename Polytope::MT MT; + typedef typename Polytope::MT MT; + typedef typename Polytope::VT VT; typedef typename Point::FT NT; - template + template Walk(GenericPolytope& P, Point const& p, - Ellipsoid const& E, // ellipsoid representing the Gaussian distribution + MT 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); + Eigen::LLT lltOfE(E.llt().solve(MT::Identity(E.cols(), E.cols()))); // compute the Cholesky decomposition of inv(E) + if (lltOfE.info() != Eigen::Success) { + throw std::runtime_error("Cholesky decomposition failed!"); + } + _L_cov = lltOfE.matrixL(); + _E = E; + _AA.noalias() = P.get_mat() * P.get_mat().transpose(); + _rho = 1000 * P.dimension(); // upper bound for the number of reflections (experimental) + initialize(P, p, rng); } - template + template Walk(GenericPolytope& P, Point const& p, - Ellipsoid const& E, // ellipsoid representing the Gaussian distribution + MT 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); + Eigen::LLT lltOfE(E.llt().solve(MT::Identity(E.cols(), E.cols()))); // compute the Cholesky decomposition of inv(E) + if (lltOfE.info() != Eigen::Success) { + throw std::runtime_error("Cholesky decomposition failed!"); + } + _L_cov = lltOfE.matrixL(); + _E = E; + _AA.noalias() = P.get_mat() * P.get_mat().transpose(); + _rho = 1000 * P.dimension(); // upper bound for the number of reflections (experimental) + initialize(P, p, rng); } - 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 +127,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 +179,42 @@ 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() = P.get_mat() * _E; + _AEA = _AE.cwiseProduct(P.get_mat()).rowwise().sum(); + + _v = GetDirection::apply(n, rng, false); + _v = Point(_L_cov.template triangularView() * _v.getCoefficients()); - NT T = -std::log(rng.sample_urdist()) * _Len; + 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 +226,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 +248,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 + MT _AA; + MT _L_cov; // LL' = inv(E) + MT _AE; + MT _E; + VT _AEA; + unsigned int _rho; update_parameters _update_parameters; typename Point::Coeff _lambdas; typename Point::Coeff _Av; + bool was_reset; }; }; diff --git a/include/random_walks/uniform_accelerated_billiard_walk.hpp b/include/random_walks/uniform_accelerated_billiard_walk.hpp index 980c8e5da..fa9c6fb8f 100644 --- a/include/random_walks/uniform_accelerated_billiard_walk.hpp +++ b/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 @@ -49,10 +50,10 @@ struct AcceleratedBilliardWalk template - < - typename Polytope, - typename RandomNumberGenerator - > + < + typename Polytope, + typename RandomNumberGenerator + > struct Walk { typedef typename Polytope::PointType Point; @@ -62,6 +63,9 @@ struct AcceleratedBilliardWalk 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); @@ -74,6 +78,9 @@ 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 @@ -104,8 +111,14 @@ struct AcceleratedBilliardWalk Point p0 = _p; it = 0; - std::pair pbpair = P.line_positive_intersect(_p, _v, _lambdas, _Av, - _lambda_prev, _update_parameters); + 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; @@ -134,7 +147,10 @@ struct AcceleratedBilliardWalk P.compute_reflection(_v, _p, _update_parameters); it++; } - if (it == _rho) _p = p0; + if (it == _rho) { + _p = p0; + was_reset = true; + } } p = _p; } @@ -185,8 +201,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 +211,7 @@ struct AcceleratedBilliardWalk } if (Lmax > _L) { - if (P.dimension() <= 2500) - { + if (P.dimension() <= 2500) { update_delta(Lmax); } else{ @@ -231,6 +245,7 @@ struct AcceleratedBilliardWalk { 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; @@ -301,6 +316,7 @@ struct AcceleratedBilliardWalk update_parameters _update_parameters; typename Point::Coeff _lambdas; typename Point::Coeff _Av; + bool was_reset; }; }; diff --git a/include/random_walks/uniform_accelerated_billiard_walk_parallel.hpp b/include/random_walks/uniform_accelerated_billiard_walk_parallel.hpp index 1c2630f9c..7ec437b23 100644 --- a/include/random_walks/uniform_accelerated_billiard_walk_parallel.hpp +++ b/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/include/sampling/random_point_generators.hpp b/include/sampling/random_point_generators.hpp index e70e14617..e60181969 100644 --- a/include/sampling/random_point_generators.hpp +++ b/include/sampling/random_point_generators.hpp @@ -95,7 +95,7 @@ struct MultivariateGaussianRandomPointGenerator Walk walk(P, p, E, rng, parameters); for (unsigned int i=0; i) add_test(NAME test_mmcs COMMAND mmcs_test -tc=mmcs) diff --git a/test/sampling_test.cpp b/test/sampling_test.cpp index 8e82e5818..887d8489d 100644 --- a/test/sampling_test.cpp +++ b/test/sampling_test.cpp @@ -28,6 +28,8 @@ #include "diagnostics/univariate_psrf.hpp" +#include "preprocess/inscribed_ellipsoid_rounding.hpp" + template < @@ -302,6 +304,53 @@ void call_test_ghmc(){ CHECK(score.maxCoeff() < 2.2); } +template +void call_test_gabw(){ + typedef Cartesian Kernel; + typedef typename Kernel::Point Point; + typedef HPolytope Hpolytope; + typedef Eigen::Matrix MT; + typedef Eigen::Matrix VT; + Hpolytope P; + unsigned int d = 10; + typedef BoostRandomNumberGenerator RNGType; + + unsigned int walkL = 10, numpoints = 10000, nburns = 0; + RNGType rng(d); + Point StartingPoint(d); + std::list randPoints; + + std::cout << "--- Testing Gaussian Accelerated Billiard Walk for Skinny-H-cube10" << std::endl; + P = generate_skinny_cube(10); + + + Point p = P.ComputeInnerBall().first; + typedef typename GaussianAcceleratedBilliardWalk::template Walk + < + Hpolytope, + RNGType + > walk; + typedef MultivariateGaussianRandomPointGenerator RandomPointGenerator; + PushBackWalkPolicy push_back_policy; + + std::tuple ellipsoid = compute_inscribed_ellipsoid + (P.get_mat(), P.get_vec(), p.getCoefficients(), 500, std::pow(10, -6.0), std::pow(10, -4.0)); + const MT E = get<0>(ellipsoid); + + RandomPointGenerator::apply(P, p, E, numpoints, 1, randPoints, + push_back_policy, rng); + + MT samples(d, numpoints); + unsigned int jj = 0; + for (typename std::list::iterator rpit = randPoints.begin(); rpit!=randPoints.end(); rpit++, jj++) + samples.col(jj) = (*rpit).getCoefficients(); + + VT score = univariate_psrf(samples); + std::cout << "psrf = " << score.maxCoeff() << std::endl; + + CHECK(score.maxCoeff() < 1.1); +} + TEST_CASE("dikin") { call_test_dikin(); } @@ -333,3 +382,7 @@ TEST_CASE("gbaw") { TEST_CASE("ghmc") { call_test_ghmc(); } + +TEST_CASE("gabw") { + call_test_gabw(); +}