Refactoring polytopes

CMakeLists.txt

@@ -0,0 +1,160 @@
+# VolEsti (volume computation and sampling library)
+# Copyright (c) 20012-2020 Vissarion Fisikopoulos
+# Copyright (c) 2018-2020 Apostolos Chalkis
+# Licensed under GNU LGPL.3, see LICENCE file
+
+project (VolEsti)
+
+enable_testing ()
+
+CMAKE_MINIMUM_REQUIRED (VERSION 2.4.5)
+
+set (CMAKE_ALLOW_LOOSE_LOOP_CONSTRUCTS true)
+
+if (COMMAND cmake_policy)
+    cmake_policy (SET CMP0003 NEW)
+endif (COMMAND cmake_policy)
+
+#if(NOT DEFINED BOOST)
+#  message(FATAL_ERROR "This program requires the boost library, and will not be compiled.")
+#else()
+#  message(STATUS "BOOST is DEFINED")
+#endif()
+
+find_library (LP_SOLVE NAMES liblpsolve55.so liblpsolve55.dylib PATHS /usr/lib/lp_solve /opt/local/lib)
+
+if (NOT LP_SOLVE)
+
+    message (FATAL_ERROR "This program requires the lp_solve library, and will not be compiled.")
+
+else ()
+
+    message (STATUS "Library lp_solve found: ${LP_SOLVE}")
+
+    set (CMAKE_EXPORT_COMPILE_COMMANDS "ON")
+
+    include_directories (BEFORE ./external/Eigen)
+    include_directories (BEFORE ./external)
+    include_directories (BEFORE ./external/minimum_ellipsoid)
+    #include_directories (BEFORE ./include/cartesian_geom)
+    #include_directories (BEFORE ./include/convex_bodies)
+    include_directories (BEFORE ./external/LPsolve_src/run_headers)
+    include_directories (BEFORE ./external/boost)
+    #include_directories (BEFORE BOOST)
+    include_directories (BEFORE ./include/generators)
+    include_directories (BEFORE ./include/volume)
+    include_directories (BEFORE ./include)
+    include_directories (BEFORE ./include/convex_bodies)
+    include_directories (BEFORE ./include/annealing)
+    include_directories (BEFORE ./include/samplers)
+    include_directories (BEFORE ./include/lp_oracles)
+    include_directories (BEFORE ./include/misc)
+
+    #for Eigen
+    if (${CMAKE_VERSION} VERSION_LESS "3.12.0")
+        add_compile_options (-D "EIGEN_NO_DEBUG")
+    else ()
+        add_compile_definitions ("EIGEN_NO_DEBUG")
+    endif ()
+
+    add_definitions (${CMAKE_CXX_FLAGS} "-std=c++11")  # enable C++11 standard
+    add_definitions (${CMAKE_CXX_FLAGS} "-g")  # enable debugger
+    #add_definitions(${CMAKE_CXX_FLAGS} "-Wint-in-bool-context")
+    # add_definitions(${CMAKE_CXX_FLAGS} "-Wall")
+
+#     add_definitions (${CMAKE_CXX_FLAGS} "-O3")  # optimization of the compiler
+    #add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lgsl")
+    add_definitions (${CXX_COVERAGE_COMPILE_FLAGS} "-lm")
+    add_definitions (${CXX_COVERAGE_COMPILE_FLAGS} "-ldl")
+    add_definitions (${CXX_COVERAGE_COMPILE_FLAGS} "-DBOOST_NO_AUTO_PTR")
+    #add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lgslcblas")
+    #add_definitions( "-O3 -lgsl -lm -ldl -lgslcblas" )
+
+
+    add_executable (generate ./test/generator.cpp)
+    TARGET_LINK_LIBRARIES (generate ${LP_SOLVE})
+    add_executable (test_generators ./test/test_generators.cpp)
+    TARGET_LINK_LIBRARIES (test_generators ${LP_SOLVE})
+
+    add_executable (new_volume_example ./test/new_volume_example.cpp)
+
+    add_executable (benchmarks_sob ./test/benchmarks_sob.cpp)
+    TARGET_LINK_LIBRARIES (benchmarks_sob ${LP_SOLVE})
+
+    add_executable (benchmarks_cg ./test/benchmarks_cg.cpp)
+    TARGET_LINK_LIBRARIES (benchmarks_cg ${LP_SOLVE})
+
+    add_executable (benchmarks_cb ./test/benchmarks_cb.cpp)
+    TARGET_LINK_LIBRARIES (benchmarks_cb ${LP_SOLVE})
+    #
+    #
+    add_library (test_main OBJECT ./test/test_main.cpp)
+
+    add_executable (volume_cg_hpolytope ./test/volume_cg_hpolytope.cpp $<TARGET_OBJECTS:test_main>)
+    add_test (NAME volume_cg_hpolytope_cube COMMAND volume_cg_hpolytope -tc=cube)
+    add_test (NAME volume_cg_hpolytope_cross COMMAND volume_cg_hpolytope -tc=cross)
+    add_test (NAME volume_cg_hpolytope_birkhoff COMMAND volume_cg_hpolytope -tc=birk)
+    add_test (NAME volume_cg_hpolytope_prod_simplex COMMAND volume_cg_hpolytope -tc=prod_simplex)
+    add_test (NAME volume_cg_hpolytope_simplex COMMAND volume_cg_hpolytope -tc=simplex)
+    add_test (NAME volume_cg_hpolytope_skinny_cube COMMAND volume_cg_hpolytope -tc=skinny_cube)
+    #
+    add_executable (volume_sob_hpolytope ./test/volume_sob_hpolytope.cpp $<TARGET_OBJECTS:test_main>)
+    add_test (NAME volume_sob_hpolytope_cube COMMAND volume_sob_hpolytope -tc=cube)
+    add_test (NAME volume_sob_hpolytope_cross COMMAND volume_sob_hpolytope -tc=cross)
+    add_test (NAME volume_sob_hpolytope_birkhoff COMMAND volume_sob_hpolytope -tc=birk)
+    add_test (NAME volume_sob_hpolytope_prod_simplex COMMAND volume_sob_hpolytope -tc=prod_simplex)
+    add_test (NAME volume_sob_hpolytope_simplex COMMAND volume_sob_hpolytope -tc=simplex)
+    add_test (NAME volume_sob_hpolytope_skinny_cube COMMAND volume_sob_hpolytope -tc=skinny_cube)
+
+
+    add_executable (volume_cb_hpolytope ./test/volume_cb_hpolytope.cpp $<TARGET_OBJECTS:test_main>)
+    add_test (NAME volume_cb_hpolytope_cube COMMAND volume_cb_hpolytope -tc=cube)
+    add_test (NAME volume_cb_hpolytope_cross COMMAND volume_cb_hpolytope -tc=cross)
+    add_test (NAME volume_cb_hpolytope_birkhoff COMMAND volume_cb_hpolytope -tc=birk)
+    add_test (NAME volume_cb_hpolytope_prod_simplex COMMAND volume_cb_hpolytope -tc=prod_simplex)
+    add_test (NAME volume_cb_hpolytope_simplex COMMAND volume_cb_hpolytope -tc=simplex)
+    add_test (NAME volume_cb_hpolytope_skinny_cube COMMAND volume_cb_hpolytope -tc=skinny_cube)
+
+    add_executable (volume_cg_vpolytope ./test/volume_cg_vpolytope.cpp $<TARGET_OBJECTS:test_main>)
+    add_test (NAME volume_cg_vpolytope_cube COMMAND volume_cg_vpolytope -tc=cube)
+    add_test (NAME volume_cg_vpolytope_cross COMMAND volume_cg_vpolytope -tc=cross)
+    add_test (NAME volume_cg_vpolytope_simplex COMMAND volume_cg_vpolytope -tc=simplex)
+
+    add_executable (volume_sob_vpolytope ./test/volume_sob_vpolytope.cpp $<TARGET_OBJECTS:test_main>)
+    add_test (NAME volume_sob_vpolytope_cube COMMAND volume_sob_vpolytope -tc=cube)
+    add_test (NAME volume_sob_vpolytope_cross COMMAND volume_sob_vpolytope -tc=cross)
+    add_test (NAME volume_sob_vpolytope_simplex COMMAND volume_sob_vpolytope -tc=simplex)
+
+
+    add_executable (volume_cb_vpolytope ./test/volume_cb_vpolytope.cpp $<TARGET_OBJECTS:test_main>)
+    add_test (NAME volume_cb_vpolytope_cube COMMAND volume_cb_vpolytope -tc=cube)
+    add_test (NAME volume_cb_vpolytope_cross COMMAND volume_cb_vpolytope -tc=cross)
+    add_test (NAME volume_cb_vpolytope_simplex COMMAND volume_cb_vpolytope -tc=simplex)
+
+    add_executable (volume_cb_zonotopes ./test/volume_cb_zonotopes.cpp $<TARGET_OBJECTS:test_main>)
+    add_test (NAME volume_cb_zonotopes_uniform_zonotopes
+              COMMAND volume_cb_zonotopes -tc=uniform_zonotopes)
+
+    add_executable (volume_cb_vpoly_intersection_vpoly ./test/volume_cb_vpoly_intersection_vpoly.cpp $<TARGET_OBJECTS:test_main>)
+    add_test (NAME volume_cb_vpoly_intersection_vpoly_random_vpoly_sphere
+              COMMAND volume_cb_vpoly_intersection_vpoly -tc=random_vpoly_sphere)
+
+    add_executable (new_rounding_test ./test/new_rounding_test.cpp $<TARGET_OBJECTS:test_main>)
+    add_test (NAME new_rounding_test_round_skinny_cube
+              COMMAND new_rounding_test -tc=round_skinny_cube)
+
+    TARGET_LINK_LIBRARIES (volume_cg_hpolytope ${LP_SOLVE})
+
+    TARGET_LINK_LIBRARIES (new_volume_example ${LP_SOLVE})
+    TARGET_LINK_LIBRARIES (volume_sob_hpolytope ${LP_SOLVE})
+    TARGET_LINK_LIBRARIES (volume_sob_vpolytope ${LP_SOLVE})
+    TARGET_LINK_LIBRARIES (volume_cg_vpolytope ${LP_SOLVE})
+    TARGET_LINK_LIBRARIES (volume_cb_hpolytope ${LP_SOLVE})
+    TARGET_LINK_LIBRARIES (volume_cb_vpolytope ${LP_SOLVE})
+    TARGET_LINK_LIBRARIES (volume_cb_zonotopes ${LP_SOLVE})
+    TARGET_LINK_LIBRARIES (volume_cb_vpoly_intersection_vpoly ${LP_SOLVE})
+    TARGET_LINK_LIBRARIES (new_rounding_test ${LP_SOLVE})
+
+endif ()
+
+

include/convex_bodies/hpolytope.h

@@ -16,38 +16,97 @@
 #include <Eigen/Eigen>
 #include "lp_oracles/solve_lp.h"
 
-// H-polytope class
+//! H-polytope class
+/*!
+    A class for a polytopne in H-representaion
+ */
 template <typename Point>
-class HPolytope{
+class HPolytope {
 public:
-    typedef Point PointType;
-    typedef typename Point::FT NT;
-    typedef typename std::vector<NT>::iterator viterator;
+    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 Eigen::Matrix<NT,Eigen::Dynamic,1> VT;
+    typedef Eigen::Matrix<NT, Eigen::Dynamic, Eigen::Dynamic> MT;
+    typedef Eigen::Matrix<NT, Eigen::Dynamic, 1>              VT;
 
 private:
-    MT A; //matrix A
-    VT b; // vector b, s.t.: Ax<=b
-    unsigned int            _d; //dimension
-    std::pair<Point,NT> _inner_ball;
-    NT maxNT = std::numeric_limits<NT>::max();
-    NT minNT = std::numeric_limits<NT>::lowest();
+    MT                   A; //matrix A
+    VT                   b; // vector b, s.t.: Ax<=b
+    unsigned int         _d; //dimension
+    std::pair<Point, NT> _inner_ball;
+    // TODO: Why the following are not static or outside the class?
+    NT                   maxNT = std::numeric_limits<NT>::max();
+    NT                   minNT = std::numeric_limits<NT>::lowest();
 
 public:
-
+    //TODO: the default implementation of the Big3 should be ok. Recheck.
     HPolytope() {}
 
-    std::pair<Point,NT> InnerBall() const
+    HPolytope(unsigned d_, const MT& A_, const VT& b_) :
+        _d{d_}, A{A_}, b{b_}, _inner_ball{ComputeChebychevBall<NT, Point>(A, b)}
+    {
+    }
+    // Copy constructor
+    HPolytope(const HPolytope<Point>& p) :
+            _d{p._d}, A{p.A}, b{p.b},  _inner_ball{p._inner_ball}
+    {
+    }
+
+    HPolytope(HPolytope&& p) :_d{p._d}
+    {
+        b = std::move(p.b);
+        _inner_ball = std::move(p._inner_ball);
+    }
+
+    HPolytope& operator=(const HPolytope& p)
+    {
+        if (this != &p) { // protect against invalid self-assignment
+            _d = p._d;
+            A = p.A;
+            b = p.b;
+            _inner_ball = p._inner_ball;
+        }
+        return *this;
+    }
+
+    HPolytope& operator=(HPolytope&& p)
+    {
+//        std::cout << __FUNCTION__ << "\n";
+        if (this != &p) {
+            _d = p._d;
+            A = std::move(p.A);
+            b = std::move(p.b);
+            _inner_ball = std::move(p._inner_ball);
+        }
+        return *this;
+    }
+
+    //define matrix A and vector b, s.t. Ax<=b and the dimension
+    HPolytope(const std::vector<std::vector<NT>>& Pin)
+    {
+        _d = Pin[0][1] - 1;
+        A.resize(Pin.size() - 1, _d);
+        b.resize(Pin.size() - 1);
+        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];
+            }
+        }
+        _inner_ball = ComputeChebychevBall<NT, Point>(A, b);
+    }
+
+
+    std::pair<Point, NT> InnerBall() const
     {
         return _inner_ball;
     }
 
     //Compute Chebyshev ball of H-polytope P:= Ax<=b
     //Use LpSolve library
-    std::pair<Point,NT> ComputeInnerBall()
+    std::pair<Point, NT> ComputeInnerBall()
     {
         _inner_ball = ComputeChebychevBall<NT, Point>(A, b);
         return _inner_ball;
@@ -109,31 +168,11 @@ class HPolytope{
         return 0.0;
     }
 
-    void init(unsigned int const& dim, MT const& _A, VT const& _b)
-    {
-        _d = dim;
-        A = _A;
-        b = _b;
-    }
-
-    //define matrix A and vector b, s.t. Ax<=b and the dimension
-    void init(std::vector<std::vector<NT> > const& Pin)
-    {
-        _d = Pin[0][1] - 1;
-        A.resize(Pin.size() - 1, _d);
-        b.resize(Pin.size() - 1);
-        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];
-            }
-        }
-    }
 
 
     // print polytope in input format
     void print() {
-        std::cout << " " << A.rows() << " " << _d << " float" << std::endl;
+        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) << " ";
@@ -257,7 +296,7 @@ class HPolytope{
             sum_nom_data++;
             sum_denom_data++;
         }
-        return std::pair<NT, NT>(min_plus, max_minus);
+        return std::make_pair(min_plus, max_minus);
     }
 
 
@@ -296,8 +335,8 @@ class HPolytope{
             Av_data++;
             sum_nom_data++;
         }
-        if (pos) return std::pair<NT, NT>(min_plus, facet);
-        return std::pair<NT, NT>(min_plus, max_minus);
+        if (pos) return std::make_pair(min_plus, facet);
+        return std::make_pair(min_plus, max_minus);
     }
 
     std::pair<NT,NT> line_intersect(Point const& r,
@@ -336,8 +375,8 @@ class HPolytope{
             Av_data++;
             sum_nom_data++;
         }
-        if (pos) return std::pair<NT, NT>(min_plus, facet);
-        return std::pair<NT, NT>(min_plus, max_minus);
+        if (pos) return std::make_pair(min_plus, facet);
+        return std::make_pair(min_plus, max_minus);
     }
 
 
@@ -396,7 +435,7 @@ class HPolytope{
             lamda_data++;
             sum_denom_data++;
         }
-        return std::pair<NT, NT>(min_plus, max_minus);
+        return std::make_pair(min_plus, max_minus);
     }
 
 
@@ -429,7 +468,7 @@ class HPolytope{
             }
             data++;
         }
-        return std::pair<NT, NT>(min_plus, max_minus);
+        return std::make_pair(min_plus, max_minus);
     }
 
 
@@ -488,8 +527,6 @@ class HPolytope{
         v += -2 * v.dot(A.row(facet)) * A.row(facet);
     }
 
-    void free_them_all() {}
-
 };
 
 #endif

include/convex_bodies/vpolyintersectvpoly.h

@@ -31,7 +31,7 @@ class IntersectionOfVpoly {
     VPolytope P1;
     VPolytope P2;
 
-    IntersectionOfVpoly() {}
+    IntersectionOfVpoly(): P1(), P2() {}
 
     IntersectionOfVpoly(VPolytope &P, VPolytope &Q) : P1(P), P2(Q) {};