From deb4cc68c41b375c65def09c0268081f2a3885a8 Mon Sep 17 00:00:00 2001 From: lucaperju <104904648+lucaperju@users.noreply.github.com> Date: Thu, 18 Jul 2024 00:28:17 +0300 Subject: [PATCH] Order polytope additions (#322) * Order Polytopes generation * minor changes for PR * remove space and comment * remove bug * Unit test for Random Order Polytope, and minor changes * remove comment --- examples/crhmc_sampling/.gitignore | 1 + include/convex_bodies/orderpolytope.h | 6 ++ include/generators/order_polytope_generator.h | 82 ++++++++++++++++++- include/misc/poset.h | 65 ++++++++------- test/order_polytope.cpp | 14 ++++ 5 files changed, 136 insertions(+), 32 deletions(-) diff --git a/examples/crhmc_sampling/.gitignore b/examples/crhmc_sampling/.gitignore index 60b43389a..1819f71c8 100644 --- a/examples/crhmc_sampling/.gitignore +++ b/examples/crhmc_sampling/.gitignore @@ -5,3 +5,4 @@ CRHMC_SIMD_* sampling_functions simple_crhmc libQD_LIB.a +liblp_solve.a diff --git a/include/convex_bodies/orderpolytope.h b/include/convex_bodies/orderpolytope.h index 79f41e161..3aa9fe8ed 100644 --- a/include/convex_bodies/orderpolytope.h +++ b/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 { diff --git a/include/generators/order_polytope_generator.h b/include/generators/order_polytope_generator.h index bfece556e..3cd907bd9 100644 --- a/include/generators/order_polytope_generator.h +++ b/include/generators/order_polytope_generator.h @@ -9,11 +9,26 @@ #ifndef ORDER_POLYTOPES_GEN_H #define ORDER_POLYTOPES_GEN_H +#include #include +#include #include -#include "misc.h" +#include +#include + +#include "misc/misc.h" #include "misc/poset.h" +#include +#include +#include +#include + +#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 instances = @@ -38,6 +53,23 @@ static const std::unordered_map instances = }; +// generates a Polytope from a poset +/// @tparam Polytope Type of returned polytope +template +Polytope get_orderpoly(Poset const &poset) { + typedef typename Polytope::PointType Point; + + OrderPolytope OP(poset); + if constexpr (std::is_same< Polytope, OrderPolytope >::value ) { + return OP; + } else if constexpr (std::is_same >::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 @@ -45,7 +77,7 @@ template 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(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(poset); +} + +// Generates a random Order Polytope with given dimension and number of facets +/// @tparam Polytope Type of returned polytope +/// @tparam RNGType RNGType Type +template +Polytope random_orderpoly(unsigned int dim, unsigned int m, int seed = std::numeric_limits::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 RNG; + RNG rng(dim); + rng.set_seed(rng_seed); + + + std::vector 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(poset); } #endif diff --git a/include/misc/poset.h b/include/misc/poset.h index 69ba02ec7..084bf66aa 100644 --- a/include/misc/poset.h +++ b/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 &res) + { + std::vector > adjList(n); + std::vector indegree(n, 0); + + for(auto x: relations) { + adjList[x.first].push_back(x.second); + indegree[x.second] += 1; + } + + std::queue q; + for(unsigned int i=0; i 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 topologically_sorted_list() const { - std::vector > adjList(n); - std::vector indegree(n, 0); - - for(auto x: order_relations) { - adjList[x.first].push_back(x.second); - indegree[x.second] += 1; - } - - std::queue q; - for(unsigned int i=0; i res; - while(!q.empty()) { - unsigned int curr = q.front(); - res.push_back(curr); - q.pop(); - - for(unsigned int i=0; i= 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 HP = random_orderpoly, NT>(10, 30); + + d = HP.dimension(); + m = HP.num_of_hyperplanes(); + + CHECK(d == 10); + CHECK(m == 30); + }