Support for sparse Eigen matrices in Hpolytope class (#327)

* 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

.gitignore

@@ -9,3 +9,4 @@ test/Testing/Temporary/CTestCostData.txt
 build/
 .vscode
 .DS_Store
+.cache

include/convex_bodies/hpolytope.h

@@ -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.
@@ -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}
     {
     }
 
@@ -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);
     }
 
@@ -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;
     }
 
@@ -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
@@ -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;
         }
@@ -493,7 +524,7 @@ class HPolytope {
                                                      VT& Ar,
                                                      VT& Av,
                                                      NT const& lambda_prev,
-                                                     MT const& AA,
+                                                     DenseMT const& AA,
                                                      update_parameters& params) const
     {
 
@@ -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;
 
@@ -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;
@@ -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;
     }
 
 
@@ -838,6 +875,7 @@ class HPolytope {
     void shift(const VT &c)
     {
         b -= A*c;
+        has_ball = false;
     }
 
 
@@ -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;
     }
@@ -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();

include/convex_bodies/orderpolytope.h

@@ -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;
     }

include/generators/h_polytopes_generator.h

@@ -10,6 +10,7 @@
 
 #include <exception>
 #include <chrono>
+#include <Eigen/Eigen>
 #include <boost/random/normal_distribution.hpp>
 #include <boost/random/uniform_real_distribution.hpp>
 
@@ -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();
@@ -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;
 

include/generators/known_polytope_generators.h

@@ -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;
@@ -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;
@@ -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;
@@ -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;
@@ -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;
@@ -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);

include/generators/order_polytope_generator.h

@@ -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";

include/lp_oracles/solve_lp.h

@@ -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);

include/random_walks/compute_diameter.hpp

@@ -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;
     }