Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Support for sparse Eigen matrices in Hpolytope class #327

Merged
merged 12 commits into from
Jul 30, 2024
Merged

Support for sparse Eigen matrices in Hpolytope class #327

Show file tree
Hide file tree
Changes from 9 commits
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
933fe04
indentation and start of sparse support
lucaperju Jul 20, 2024
29ca44e
sparse matrix support for hpolytope and gabw
lucaperju Jul 22, 2024
83a38be
unit tests and minor fixes
lucaperju Jul 22, 2024
2dea296
fix conflict bug
lucaperju Jul 22, 2024
2ad4bfa
remove debug print and change variable names
lucaperju Jul 22, 2024
4498490
fix innerBall bug
lucaperju Jul 23, 2024
bf93e1f
has_ball flag, separate A_type and E_type in GABW, improve unit test
lucaperju Jul 24, 2024
08c0f66
rounding and volume_cg unit tests for sparse hpoly
lucaperju Jul 27, 2024
fe8f92c
sampling test for sparse polytope
lucaperju Jul 27, 2024
0394b98
minor changes
lucaperju Jul 30, 2024
6fe2cf5
fix bug, update set_interior_point function
lucaperju Jul 30, 2024
52c731e
check that interior point is valid
lucaperju Jul 30, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions .gitignore
View file
Open in desktop
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
105 changes: 65 additions & 40 deletions include/convex_bodies/hpolytope.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -41,23 +41,29 @@ 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;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do we need those two lines?

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 +74,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 +96,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,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)
Copy link
Member

@TolisChal TolisChal Jul 29, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you add a TODO here to fix this function?
Because when we change the center of the inner ball we also need to change its radius.
Also the has_ball should be set to false here, right?

Copy link
Contributor Author

@lucaperju lucaperju Jul 30, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I decided to set the has_ball to true here since I assumed that whenever this function might be used, it would be given a valid innerball, but tbh I haven't seen anywhere in the code where this function is called anyways.
edit: oh, if you're refering the the set_interior_point function, then yeah, has_ball should likely be set to false, but I don't really know how we could fix this function further than that

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fisrt, we need to check that the point is indeed inside the polytope and only then set it as center of the inner ball. Then, we need to compute the distances of the center from all the facets and take the minimum. The latter distance is the radius of the maximum inscribed ball centered at the given point.

The distance between a point p and the i-th facet is given by: (b_i - a_i^Tp) / ||a_i||.
Check also the function get_dists() below; it returns the distance from the origin to each facet.

Expand All @@ -112,29 +125,29 @@ class HPolytope {
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 +198,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 +225,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 @@ -509,7 +524,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() += (DenseMT)((-2.0 * inner_prev) * AA.col(params.facet_prev));
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Input matrix AA should always be dense. Then, (DenseMT) is not needed here.

}
sum_nom.noalias() = b - Ar;

Expand Down Expand Up @@ -630,12 +645,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)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Would lamdas.noalias() += (VT)(A.col(rand_coord_prev) ... work here as well?

* (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 +841,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 +859,7 @@ class HPolytope {
void shift(const VT &c)
{
b -= A*c;
has_ball = false;
}


Expand Down Expand Up @@ -867,12 +889,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 +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
4 changes: 2 additions & 2 deletions include/convex_bodies/orderpolytope.h
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Loading