Skip to content

Commit

Permalink
Support for sparse Eigen matrices in Hpolytope class (#327)
Browse files Browse the repository at this point in the history 
* indentation and start of sparse support

* sparse matrix support for hpolytope and gabw

* unit tests and minor fixes

* fix conflict bug

* remove debug print and change variable names

* fix innerBall bug

* has_ball flag, separate A_type and E_type in GABW, improve unit test

* rounding and volume_cg unit tests for sparse hpoly

* sampling test for sparse polytope

* minor changes

* fix bug, update set_interior_point function

* check that interior point is valid
  • Loading branch information
@lucaperju
lucaperju authored Jul 30, 2024
1 parent ad5b3bf commit d424f3e
Show file tree
Hide file tree
Showing 14 changed files with 326 additions and 92 deletions.
1 change: 1 addition & 0 deletions .gitignore
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,3 +9,4 @@ test/Testing/Temporary/CTestCostData.txt
build/
.vscode
.DS_Store
.cache
129 changes: 85 additions & 44 deletions include/convex_bodies/hpolytope.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,23 +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 <typename Point>
template
<
typename Point,
typename MT_type = Eigen::Matrix<typename Point::FT, Eigen::Dynamic, Eigen::Dynamic>
>
class HPolytope {
public:
typedef Point PointType;
typedef typename Point::FT NT;
typedef typename std::vector<NT>::iterator viterator;
//using RowMatrixXd = Eigen::Matrix<NT, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
//typedef RowMatrixXd MT;
typedef Eigen::Matrix<NT, Eigen::Dynamic, Eigen::Dynamic> MT;
typedef MT_type MT;
typedef Eigen::Matrix<NT, Eigen::Dynamic, 1> VT;
typedef Eigen::Matrix<NT, Eigen::Dynamic, Eigen::Dynamic> DenseMT;

private:
unsigned int _d; //dimension
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
bool has_ball = false;

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

template<typename T = DenseMT>
HPolytope(unsigned d_, DenseMT const& A_, VT const& b_, typename std::enable_if<!std::is_same<MT, T>::value, T>::type* = 0) :
_d{d_}, A{A_.sparseView()}, b{b_}
{
}

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

Expand All @@ -84,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<NT, Point>(A, b);
}

Expand All @@ -100,41 +110,60 @@ 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)
{
_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<NT>::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<Point, NT> ComputeInnerBall()
{
normalize();
#ifndef DISABLE_LPSOLVE
_inner_ball = ComputeChebychevBall<NT, Point>(A, b); // use lpsolve library
#else

if (_inner_ball.second <= NT(0)) {

NT const tol = 1e-08;
std::tuple<VT, NT, bool> 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<VT, NT, bool> 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<NT, Point>(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;
}

Expand Down Expand Up @@ -185,13 +214,15 @@ class HPolytope {
{
A = A2;
normalized = false;
has_ball = false;
}


// change the vector b
void set_vec(VT const& b2)
{
b = b2;
has_ball = false;
}

Point get_mean_of_vertices() const
Expand All @@ -210,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;
}
Expand Down Expand Up @@ -493,7 +524,7 @@ class HPolytope {
VT& Ar,
VT& Av,
NT const& lambda_prev,
MT const& AA,
DenseMT const& AA,
update_parameters& params) const
{

Expand All @@ -509,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;

Expand Down Expand Up @@ -630,12 +661,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"<<std::endl;
Expand Down Expand Up @@ -826,10 +857,16 @@ class HPolytope {


// Apply linear transformation, of square matrix T^{-1}, in H-polytope P:= Ax<=b
void linear_transformIt(MT const& T)
template<typename T_type>
void linear_transformIt(T_type const& T)
{
A = A * T;
if constexpr (std::is_same<MT, DenseMT>::value) {
A = A * T;
} else {
A = (A * T).sparseView();
}
normalized = false;
has_ball = false;
}


Expand All @@ -838,6 +875,7 @@ class HPolytope {
void shift(const VT &c)
{
b -= A*c;
has_ball = false;
}


Expand Down Expand Up @@ -867,12 +905,15 @@ 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;
}
Expand Down Expand Up @@ -919,7 +960,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
4 changes: 2 additions & 2 deletions include/convex_bodies/orderpolytope.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -94,8 +94,8 @@ class OrderPolytope {
return _A.sparseView();
}

// return the matrix A
MT get_full_mat() const
// return the matrix A
MT get_dense_mat() const
{
return _A;
}
Expand Down
5 changes: 3 additions & 2 deletions include/generators/h_polytopes_generator.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@

#include <exception>
#include <chrono>
#include <Eigen/Eigen>
#include <boost/random/normal_distribution.hpp>
#include <boost/random/uniform_real_distribution.hpp>

Expand All @@ -25,9 +26,9 @@
template <class Polytope, class RNGType>
Polytope random_hpoly(unsigned int dim, unsigned int m, int seed = std::numeric_limits<int>::signaling_NaN()) {

typedef typename Polytope::MT MT;
typedef typename Polytope::VT VT;
typedef typename Polytope::NT NT;
typedef typename Eigen::Matrix<NT, Eigen::Dynamic, Eigen::Dynamic> MT;
typedef typename Polytope::PointType Point;

int rng_seed = std::chrono::system_clock::now().time_since_epoch().count();
Expand Down Expand Up @@ -104,7 +105,7 @@ template <class Polytope, typename NT, class RNGType>
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<int>::signaling_NaN()) {

typedef typename Polytope::MT MT;
typedef typename Eigen::Matrix<NT, Eigen::Dynamic, Eigen::Dynamic> MT;
typedef typename Polytope::VT VT;
typedef typename Polytope::PointType Point;

Expand Down
12 changes: 6 additions & 6 deletions include/generators/known_polytope_generators.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@
template <class Polytope>
Polytope generate_cube(const unsigned int& dim, const bool& Vpoly,typename Polytope::NT scale=1) {

typedef typename Polytope::MT MT;
typedef typename Eigen::Matrix<typename Polytope::NT, Eigen::Dynamic, Eigen::Dynamic> MT;
typedef typename Polytope::VT VT;
MT A;
VT b;
Expand Down Expand Up @@ -84,7 +84,7 @@ template <typename Polytope>
Polytope generate_cross(const unsigned int &dim, const bool &Vpoly) {

unsigned int m;
typedef typename Polytope::MT MT;
typedef typename Eigen::Matrix<typename Polytope::NT, Eigen::Dynamic, Eigen::Dynamic> MT;
typedef typename Polytope::VT VT;

MT A;
Expand Down Expand Up @@ -145,7 +145,7 @@ Polytope generate_cross(const unsigned int &dim, const bool &Vpoly) {
/// @tparam Polytope Type of returned polytope
template <typename Polytope>
Polytope generate_simplex(const unsigned int &dim, const bool &Vpoly){
typedef typename Polytope::MT MT;
typedef typename Eigen::Matrix<typename Polytope::NT, Eigen::Dynamic, Eigen::Dynamic> MT;
typedef typename Polytope::VT VT;

MT A;
Expand Down Expand Up @@ -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<typename Polytope::NT, Eigen::Dynamic, Eigen::Dynamic> MT;
typedef typename Polytope::VT VT;

MT A;
Expand Down Expand Up @@ -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<typename Polytope::NT, Eigen::Dynamic, Eigen::Dynamic> MT;
typedef typename Polytope::VT VT;

MT A;
Expand Down Expand Up @@ -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<typename Polytope::NT, Eigen::Dynamic, Eigen::Dynamic> MT;
typedef typename Polytope::VT VT;

MT A = MT::Zero(m, d);
Expand Down
2 changes: 1 addition & 1 deletion include/generators/order_polytope_generator.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,7 @@ Polytope get_orderpoly(Poset const &poset) {
if constexpr (std::is_same< Polytope, OrderPolytope<Point> >::value ) {
return OP;
} else if constexpr (std::is_same<Polytope, HPolytope<Point> >::value ){
Polytope HP(OP.dimension(), OP.get_full_mat(), OP.get_vec());
Polytope HP(OP.dimension(), OP.get_dense_mat(), OP.get_vec());
return HP;
} else {
throw "Unable to generate an Order Polytope of requested type";
Expand Down
4 changes: 2 additions & 2 deletions include/lp_oracles/solve_lp.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -80,8 +80,8 @@ std::pair<Point,NT> ComputeChebychevBall(const MT &A, const VT &b){
sum=NT(0);
for(j=0; j<d; j++){
colno[j] = j+1;
row[j] = A(i,j);
sum+=A(i,j)*A(i,j);
row[j] = A.coeff(i,j);
sum+=A.coeff(i,j)*A.coeff(i,j);
}
colno[d] = d+1; /* last column */
row[d] = std::sqrt(sum);
Expand Down
6 changes: 3 additions & 3 deletions include/random_walks/compute_diameter.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,11 +30,11 @@ struct compute_diameter
};


template <typename Point>
struct compute_diameter<HPolytope<Point>>
template <typename Point, typename MT>
struct compute_diameter<HPolytope<Point, MT>>
{
template <typename NT>
static NT compute(HPolytope<Point> &P)
static NT compute(HPolytope<Point, MT> &P)
{
return NT(2) * std::sqrt(NT(P.dimension())) * P.InnerBall().second;
}
Expand Down
Loading

0 comments on commit d424f3e

Please sign in to comment.