Skip to content

Commit

Permalink
has_ball flag, separate A_type and E_type in GABW, improve unit test
Browse files Browse the repository at this point in the history
  • Loading branch information
@lucaperju
lucaperju committed Jul 24, 2024
1 parent 4498490 commit bf93e1f
Show file tree
Hide file tree
Showing 3 changed files with 47 additions and 39 deletions.
27 changes: 16 additions & 11 deletions include/convex_bodies/hpolytope.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,6 +63,7 @@ class HPolytope {
VT b; // vector b, s.t.: Ax<=b
std::pair<Point, NT> _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.
Expand All @@ -81,7 +82,7 @@ class HPolytope {

// Copy constructor
HPolytope(HPolytope<Point, MT> const& p) :
_d{p._d}, A{p.A}, b{p.b}, _inner_ball{p._inner_ball}, normalized{p.normalized}
_d{p._d}, A{p.A}, b{p.b}, _inner_ball{p._inner_ball}, normalized{p.normalized}, has_ball{p.has_ball}
{
}

Expand All @@ -98,7 +99,7 @@ class HPolytope {
A(i - 1, j - 1) = -Pin[i][j];
}
}
_inner_ball.second = -1;
has_ball = false;
//_inner_ball = ComputeChebychevBall<NT, Point>(A, b);
}

Expand All @@ -111,6 +112,7 @@ class HPolytope {
void set_InnerBall(std::pair<Point,NT> const& innerball) //const
{
_inner_ball = innerball;
has_ball = true;
}

void set_interior_point(Point const& r)
Expand All @@ -120,11 +122,12 @@ class HPolytope {

//Compute Chebyshev ball of H-polytope P:= Ax<=b
//Use LpSolve library
std::pair<Point, NT> ComputeInnerBall(bool const &force = false)
std::pair<Point, NT> ComputeInnerBall()
{
normalize();
if (force || _inner_ball.second <= NT(0)) {

if (!has_ball) {

has_ball = true;
NT const tol = 1e-08;
std::tuple<VT, NT, bool> inner_ball = max_inscribed_ball(A, b, 5000, tol);

Expand All @@ -138,7 +141,7 @@ class HPolytope {
_inner_ball = ComputeChebychevBall<NT, Point>(A, b); // use lpsolve library
#else
std::cerr << "lpsolve is disabled, unable to compute inner ball";
_inner_ball.second = -1.0;
has_ball = false;
#endif
} else {
_inner_ball.first = Point(std::get<0>(inner_ball));
Expand Down Expand Up @@ -195,15 +198,15 @@ class HPolytope {
{
A = A2;
normalized = false;
_inner_ball.second = -1.0;
has_ball = false;
}


// change the vector b
void set_vec(VT const& b2)
{
b = b2;
_inner_ball.second = -1.0;
has_ball = false;
}

Point get_mean_of_vertices() const
Expand Down Expand Up @@ -847,7 +850,7 @@ class HPolytope {
A = (A * T).sparseView();
}
normalized = false;
_inner_ball.second = -1.0;
has_ball = false;
}


Expand All @@ -856,7 +859,7 @@ class HPolytope {
void shift(const VT &c)
{
b -= A*c;
_inner_ball.second = -1.0;
has_ball = false;
}


Expand Down Expand Up @@ -886,6 +889,8 @@ class HPolytope {

void normalize()
{
if(normalized)
return;
NT row_norm;
for (int i = 0; i < A.rows(); ++i) {
row_norm = A.row(i).norm();
Expand Down Expand Up @@ -939,7 +944,7 @@ class HPolytope {
}

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 {
NT compute_reflection(Point &v, const Point &, DenseMT 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();
Expand Down
44 changes: 21 additions & 23 deletions include/random_walks/gaussian_accelerated_billiard_walk.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -61,34 +61,36 @@ struct GaussianAcceleratedBilliardWalk
template
<
typename Polytope,
typename RandomNumberGenerator
typename RandomNumberGenerator,
typename E_type = typename Polytope::DenseMT
>
struct Walk
{
typedef typename Polytope::PointType Point;
typedef typename Polytope::MT MT;
typedef typename Polytope::MT A_type;
typedef typename Polytope::DenseMT DenseMT;
typedef typename Polytope::VT VT;
typedef typename Point::FT NT;

void computeLcov(const MT E)
void computeLcov(const E_type E)
{
if constexpr (std::is_same<MT, Eigen::SparseMatrix<NT> >::value) {
Eigen::SimplicialLLT<MT> lltofE;
if constexpr (std::is_same<E_type, Eigen::SparseMatrix<NT> >::value) {
Eigen::SimplicialLLT<E_type> lltofE;
lltofE.compute(E);
if (lltofE.info() != Eigen::Success) {
throw std::runtime_error("First Cholesky decomposition failed for sparse matrix!");
}
Eigen::SparseMatrix<NT> I(E.cols(), E.cols());
I.setIdentity();
Eigen::SparseMatrix<NT> E_inv = lltofE.solve(I);
Eigen::SimplicialLLT<MT> lltofEinv;
Eigen::SimplicialLLT<E_type> 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<MT> lltOfE(E.llt().solve(MT::Identity(E.cols(), E.cols()))); // compute the Cholesky decomposition of inv(E)
Eigen::LLT<E_type> 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!");
}
Expand All @@ -99,7 +101,7 @@ struct GaussianAcceleratedBilliardWalk
template <typename GenericPolytope>
Walk(GenericPolytope& P,
Point const& p,
MT const& E, // covariance matrix representing the Gaussian distribution
E_type const& E, // covariance matrix representing the Gaussian distribution
RandomNumberGenerator &rng)
{
if(!P.is_normalized()) {
Expand All @@ -109,7 +111,7 @@ struct GaussianAcceleratedBilliardWalk
_L = compute_diameter<GenericPolytope>::template compute<NT>(P);
computeLcov(E);
_E = E;
if constexpr (std::is_same<MT, Eigen::SparseMatrix<NT>>::value) {
if constexpr (std::is_same<A_type, Eigen::SparseMatrix<NT>>::value) {
_AA = P.get_mat() * P.get_mat().transpose();
} else {
_AA.noalias() = P.get_mat() * P.get_mat().transpose();
Expand All @@ -121,7 +123,7 @@ struct GaussianAcceleratedBilliardWalk
template <typename GenericPolytope>
Walk(GenericPolytope& P,
Point const& p,
MT const& E, // covariance matrix representing the Gaussian distribution
E_type const& E, // covariance matrix representing the Gaussian distribution
RandomNumberGenerator &rng,
parameters const& params)
{
Expand All @@ -134,7 +136,7 @@ struct GaussianAcceleratedBilliardWalk
::template compute<NT>(P);
computeLcov(E);
_E = E;
if constexpr (std::is_same<MT, Eigen::SparseMatrix<NT>>::value) {
if constexpr (std::is_same<A_type, Eigen::SparseMatrix<NT>>::value) {
_AA = P.get_mat() * P.get_mat().transpose();
} else {
_AA.noalias() = P.get_mat() * P.get_mat().transpose();
Expand Down Expand Up @@ -230,18 +232,14 @@ struct GaussianAcceleratedBilliardWalk
_lambdas.setZero(P.num_of_hyperplanes());
_Av.setZero(P.num_of_hyperplanes());
_p = p;
if constexpr (std::is_same<MT, Eigen::SparseMatrix<NT>>::value) {
_AE = P.get_mat() * _E;
} else {
_AE.noalias() = P.get_mat() * _E;
}
//_AEA = _AE.cwiseProduct(P.get_mat()).rowwise().sum();

_AE.noalias() = (DenseMT)(P.get_mat() * _E);
_AEA = _AE.cwiseProduct((DenseMT)P.get_mat()).rowwise().sum();
/*
_AEA.resize(P.num_of_hyperplanes());
for(int i = 0; i < P.num_of_hyperplanes(); ++i)
{
_AEA(i) = _AE.row(i).dot(P.get_mat().row(i));
}
}*/

_v = GetDirection<Point>::apply(n, rng, false);
_v = Point(_L_cov.template triangularView<Eigen::Lower>() * _v.getCoefficients());
Expand Down Expand Up @@ -293,10 +291,10 @@ struct GaussianAcceleratedBilliardWalk
Point _p;
Point _v;
NT _lambda_prev;
MT _AA;
MT _L_cov; // LL' = inv(E)
MT _AE;
MT _E;
A_type _AA;
E_type _L_cov; // LL' = inv(E)
DenseMT _AE;
E_type _E;
VT _AEA;
unsigned int _rho;
update_parameters _update_parameters;
Expand Down
15 changes: 10 additions & 5 deletions test/sampling_test.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -349,7 +349,7 @@ void call_test_gabw(){
std::cout << "psrf = " << score.maxCoeff() << std::endl;

CHECK(score.maxCoeff() < 1.1);

RNGType Srng(d);

typedef Eigen::SparseMatrix<NT> SparseMT;
typedef HPolytope<Point, SparseMT> SparseHpoly;
Expand All @@ -366,7 +366,8 @@ void call_test_gabw(){
typedef typename GaussianAcceleratedBilliardWalk::template Walk
<
SparseHpoly,
RNGType
RNGType,
SparseMT
> sparsewalk;
typedef MultivariateGaussianRandomPointGenerator <sparsewalk> SparseRandomPointGenerator;

Expand All @@ -377,16 +378,20 @@ void call_test_gabw(){
const SparseMT SE = get<0>(ellipsoid).sparseView();

SparseRandomPointGenerator::apply(SP, p, SE, numpoints, 1, Points,
push_back_policy, rng);
push_back_policy, Srng);

jj = 0;
MT Ssamples(d, numpoints);
for (typename std::list<Point>::iterator rpit = Points.begin(); rpit != Points.end(); rpit++, jj++)
samples.col(jj) = (*rpit).getCoefficients();
Ssamples.col(jj) = (*rpit).getCoefficients();

score = univariate_psrf<NT, VT>(samples);
score = univariate_psrf<NT, VT>(Ssamples);
std::cout << "psrf = " << score.maxCoeff() << std::endl;

CHECK(score.maxCoeff() < 1.1);
NT delta = (samples - Ssamples).maxCoeff();
std::cout << "\ndelta = " << delta << std::endl;
CHECK(delta < 0.00001);

}

Expand Down

0 comments on commit bf93e1f

Please sign in to comment.