Skip to content

Commit

Permalink
Gaussian Accelerated Billiard Walk (#323)
Browse files Browse the repository at this point in the history
* 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
  • Loading branch information
lucaperju authored Jul 22, 2024
1 parent c2504a5 commit 9d2eb4c
Show file tree
Hide file tree
Showing 16 changed files with 274 additions and 115 deletions.
16 changes: 6 additions & 10 deletions examples/sampling-hpolytope-with-billiard-walks/sampler.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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"

Expand Down Expand Up @@ -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<Point> 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<MT>(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<MT, VT, NT> ellipsoid = compute_inscribed_ellipsoid<MT, EllipsoidType::VOLUMETRIC_BARRIER>
(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);
}
Expand Down
6 changes: 6 additions & 0 deletions include/convex_bodies/ballintersectconvex.h
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -261,6 +261,12 @@ class CorrelationSpectrahedron : public Spectrahedron<Point>{
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
25 changes: 23 additions & 2 deletions include/convex_bodies/hpolytope.h
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -56,6 +57,7 @@ class HPolytope {
MT A; //matrix A
VT b; // vector b, s.t.: Ax<=b
std::pair<Point, NT> _inner_ball;
bool normalized = false; // true if the polytope is normalized

public:
//TODO: the default implementation of the Big3 should be ok. Recheck.
Expand All @@ -68,7 +70,7 @@ class HPolytope {

// Copy constructor
HPolytope(HPolytope<Point> 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}
{
}

Expand Down Expand Up @@ -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;
}


Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -822,6 +829,7 @@ class HPolytope {
void linear_transformIt(MT const& T)
{
A = A * T;
normalized = false;
}


Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -909,6 +918,18 @@ class HPolytope {
v += a;
}

template <typename update_parameters>
NT compute_reflection(Point &v, const Point &, MT const &AE, VT const &AEA, NT const &vEv, update_parameters const &params) 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 <class bfunc, class NonLinearOracle>
Expand Down
4 changes: 4 additions & 0 deletions include/convex_bodies/orderpolytope.h
Original file line number Diff line number Diff line change
Expand Up @@ -100,6 +100,10 @@ class OrderPolytope {
return b;
}

bool is_normalized ()
{
return _normalized;
}

// print polytope in Ax <= b format
void print() const
Expand Down
7 changes: 7 additions & 0 deletions include/convex_bodies/spectrahedra/spectrahedron.h
Original file line number Diff line number Diff line change
Expand Up @@ -341,6 +341,13 @@ class Spectrahedron {
{
return MT::Identity(lmi.dimension(), lmi.dimension());
}

bool is_normalized ()
{
return true;
}

void normalize() {}

void resetFlags()
{
Expand Down
4 changes: 4 additions & 0 deletions include/convex_bodies/vpolyintersectvpoly.h
Original file line number Diff line number Diff line change
Expand Up @@ -92,6 +92,10 @@ class IntersectionOfVpoly {
return P1.get_vec();
}

bool is_normalized () {
return true;
}

MT get_T() const {
return P1.get_mat();
}
Expand Down
4 changes: 4 additions & 0 deletions include/convex_bodies/vpolytope.h
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
6 changes: 6 additions & 0 deletions include/convex_bodies/zonoIntersecthpoly.h
Original file line number Diff line number Diff line change
Expand Up @@ -65,6 +65,12 @@ class ZonoIntersectHPoly {
return HP.get_vec();
}

bool is_normalized () {
return true;
}

void normalize() {}

std::pair<PointType,NT> InnerBall() const
{
return Z.InnerBall();
Expand Down
5 changes: 5 additions & 0 deletions include/convex_bodies/zpolytope.h
Original file line number Diff line number Diff line change
Expand Up @@ -299,6 +299,11 @@ class Zonotope {
return b;
}

bool is_normalized ()
{
return true;
}

// change the matrix V
void set_mat(MT const& V2)
{
Expand Down
Loading

0 comments on commit 9d2eb4c

Please sign in to comment.