Order polytope improvements (#319)

* Order Polytopes generation

* minor changes for PR

* remove space and comment

* remove bug

* Unit test for Random Order Polytope, and minor changes

* remove comment

* Order Polytopes generation

* minor changes for PR

* remove space and comment

* remove bug

* Unit test for Random Order Polytope, and minor changes

* remove comment

* fix rebase bugs

* remove accidental line

Author: lucaperju

examples/crhmc_sampling/.gitignore

@@ -5,3 +5,4 @@ CRHMC_SIMD_*
 sampling_functions
 simple_crhmc
 libQD_LIB.a
+liblp_solve.a

include/convex_bodies/orderpolytope.h

@@ -94,6 +94,12 @@ class OrderPolytope {
         return _A.sparseView();
     }
 
+    // return the matrix A
+    MT get_full_mat() const
+    {
+        return _A;
+    }
+
 
     VT get_vec() const
     {

include/generators/order_polytope_generator.h

@@ -9,11 +9,26 @@
 #ifndef ORDER_POLYTOPES_GEN_H
 #define ORDER_POLYTOPES_GEN_H
 
+#include <chrono>
 #include <sstream>
+#include <type_traits>
 #include <unordered_map>
-#include "misc.h"
+#include <vector>
+#include <algorithm>
+
+#include "misc/misc.h"
 #include "misc/poset.h"
 
+#include <boost/random.hpp>
+#include <boost/random/uniform_int.hpp>
+#include <boost/random/normal_distribution.hpp>
+#include <boost/random/uniform_real_distribution.hpp>
+
+#include "generators/boost_random_number_generator.hpp"
+
+#include "convex_bodies/orderpolytope.h"
+#include "convex_bodies/hpolytope.h"
+
 
 // Instances taken from: https://github.com/ttalvitie/le-counting-practice
 static const std::unordered_map<std::string, std::string> instances =
@@ -38,14 +53,31 @@ static const std::unordered_map<std::string, std::string> instances =
 
 };
 
+// generates a Polytope from a poset
+/// @tparam Polytope Type of returned polytope
+template <class Polytope>
+Polytope get_orderpoly(Poset const &poset) {
+    typedef typename Polytope::PointType Point;
+
+    OrderPolytope<Point> OP(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());
+        return HP;
+    } else {
+        throw "Unable to generate an Order Polytope of requested type";
+    }
+}
+
 // generates an Order Polytope from an instance name
 // Instances taken from: https://github.com/ttalvitie/le-counting-practice
 /// @tparam Polytope Type of returned polytope
 template <class Polytope>
 Polytope generate_orderpoly(std::string& instance_name) {
     std::stringstream in_ss(instances.at(instance_name));
     Poset poset = read_poset_from_file_adj_matrix(in_ss).second;
-    return Polytope(poset);
+    return get_orderpoly<Polytope>(poset);
 }
 
 // Generates a cube as an Order Polytope
@@ -56,8 +88,50 @@ Polytope generate_cube_orderpoly(unsigned int dim) {
 
     RV order_relations;
     Poset poset(dim, order_relations);
-    Polytope OP(poset);
-    return OP;
+    return get_orderpoly<Polytope>(poset);
+}
+
+// Generates a random Order Polytope with given dimension and number of facets
+/// @tparam Polytope Type of returned polytope
+/// @tparam RNGType RNGType Type
+template <class Polytope, typename NT>
+Polytope random_orderpoly(unsigned int dim, unsigned int m, int seed = std::numeric_limits<int>::signaling_NaN()) {
+
+    typedef typename Poset::RV RV;
+
+    int rng_seed = std::chrono::system_clock::now().time_since_epoch().count();
+    if (!isnan(seed)) {
+        rng_seed = seed;
+    }
+
+    typedef BoostRandomNumberGenerator<boost::mt19937, NT> RNG;
+    RNG rng(dim);
+    rng.set_seed(rng_seed);
+
+
+    std::vector<unsigned int> order(dim);
+    for(int i = 0; i < dim; ++i) {
+        order[i] = i;
+    }
+    boost::mt19937 shuffle_rng(rng_seed);
+    std::shuffle(order.begin(), order.end(), shuffle_rng);
+
+
+    RV order_relations;
+    for(int i = 0; i < m - 2 * dim; ++i) {
+        unsigned int x = rng.sample_uidist();
+        unsigned int y = rng.sample_uidist();
+        while(x == y) {
+            y = rng.sample_uidist();
+        }
+        if(x > y)
+            std::swap(x, y);
+        order_relations.push_back(std::make_pair(order[x], order[y]));
+    }
+
+
+    Poset poset(dim, order_relations);
+    return get_orderpoly<Polytope>(poset);
 }
 
 #endif

include/misc/poset.h

@@ -25,6 +25,36 @@ class Poset {
     unsigned int n;  // elements will be from 0 to n-1
     RV order_relations;     // pairs of form a <= b
 
+    static void sorted_list(const unsigned int &n, const RV &relations, std::vector<unsigned int> &res)
+    {
+        std::vector<std::vector<unsigned int> > adjList(n);
+        std::vector<unsigned int> indegree(n, 0);
+
+        for(auto x: relations) {
+            adjList[x.first].push_back(x.second);
+            indegree[x.second] += 1;
+        }
+
+        std::queue<unsigned int> q;
+        for(unsigned int i=0; i<n; ++i) {
+            if(indegree[i] == 0)
+                q.push(i);
+        }
+
+        while(!q.empty()) {
+            unsigned int curr = q.front();
+            res.push_back(curr);
+            q.pop();
+
+            for(unsigned int i=0; i<adjList[curr].size(); ++i) {
+                unsigned int adj_idx = adjList[curr][i];
+                indegree[adj_idx] -= 1;
+                if(indegree[adj_idx] == 0)
+                    q.push(adj_idx);
+            }
+        }
+    }
+
 public:
     Poset() {}
 
@@ -44,7 +74,12 @@ class Poset {
                 throw "invalid elements in order relations";
         }
 
-        // TODO: Check if corresponding DAG is actually acyclic
+        std::vector<unsigned int> order;
+        sorted_list(n, relations, order);
+        
+        if(order.size() < n) { // TODO: accept cycles in the poset
+            throw "corresponding DAG is not acyclic";
+        }
 
         return relations;
     }
@@ -96,34 +131,8 @@ class Poset {
 
     std::vector<unsigned int> topologically_sorted_list() const
     {
-        std::vector<std::vector<unsigned int> > adjList(n);
-        std::vector<unsigned int> indegree(n, 0);
-
-        for(auto x: order_relations) {
-            adjList[x.first].push_back(x.second);
-            indegree[x.second] += 1;
-        }
-
-        std::queue<unsigned int> q;
-        for(unsigned int i=0; i<n; ++i) {
-            if(indegree[i] == 0)
-                q.push(i);
-        }
-
         std::vector<unsigned int> res;
-        while(!q.empty()) {
-            unsigned int curr = q.front();
-            res.push_back(curr);
-            q.pop();
-
-            for(unsigned int i=0; i<adjList[curr].size(); ++i) {
-                unsigned int adj_idx = adjList[curr][i];
-                indegree[adj_idx] -= 1;
-                if(indegree[adj_idx] == 0)
-                    q.push(adj_idx);
-            }
-        }
-
+        sorted_list(n, order_relations, res);
         return res;
     }
 };

test/order_polytope.cpp

@@ -21,6 +21,10 @@
 #include "cartesian_geom/cartesian_kernel.h"
 #include "cartesian_geom/point.h"
 #include "convex_bodies/orderpolytope.h"
+#include "convex_bodies/hpolytope.h"
+
+#include "generators/order_polytope_generator.h"
+
 #include "misc/poset.h"
 #include "misc/misc.h"
 
@@ -150,6 +154,16 @@ void call_test_basics() {
     CHECK(OP.is_in(Point(4, {0.5, 0.5, 0.0, 1.0})) == 0);   // a0 <= a2 violated
     CHECK(OP.is_in(Point(4, {-0.1, 0.5, 1.0, 1.0})) == 0);  // a0 >= 0 violated
     CHECK(OP.is_in(Point(4, {1.0, 0.5, 1.0, 1.1})) == 0);   // a3 <= 1 violated
+
+    // Create a random Order Polytope of dimension 10 with 30 facets as an Hpolytope class
+    HPolytope<Point> HP = random_orderpoly<HPolytope<Point>, NT>(10, 30);
+
+    d = HP.dimension();
+    m = HP.num_of_hyperplanes();
+
+    CHECK(d == 10);
+    CHECK(m == 30);
+    
 }