Skip to content

Commit

Permalink
Order polytope additions (#322)
Browse files Browse the repository at this point in the history
* Order Polytopes generation

* minor changes for PR

* remove space and comment

* remove bug

* Unit test for Random Order Polytope, and minor changes

* remove comment
  • Loading branch information
lucaperju authored Jul 17, 2024
1 parent 7a0be2d commit deb4cc6
Show file tree
Hide file tree
Showing 5 changed files with 136 additions and 32 deletions.
1 change: 1 addition & 0 deletions examples/crhmc_sampling/.gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -5,3 +5,4 @@ CRHMC_SIMD_*
sampling_functions
simple_crhmc
libQD_LIB.a
liblp_solve.a
6 changes: 6 additions & 0 deletions include/convex_bodies/orderpolytope.h
Original file line number Diff line number Diff line change
Expand Up @@ -94,6 +94,12 @@ class OrderPolytope {
return _A.sparseView();
}

// return the matrix A
MT get_full_mat() const
{
return _A;
}


VT get_vec() const
{
Expand Down
82 changes: 78 additions & 4 deletions include/generators/order_polytope_generator.h
Original file line number Diff line number Diff line change
Expand Up @@ -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 =
Expand All @@ -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
Expand All @@ -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
65 changes: 37 additions & 28 deletions include/misc/poset.h
Original file line number Diff line number Diff line change
Expand Up @@ -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() {}

Expand All @@ -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;
}
Expand Down Expand Up @@ -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;
}
};
Expand Down
14 changes: 14 additions & 0 deletions test/order_polytope.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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"

Expand Down Expand Up @@ -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);

}


Expand Down

0 comments on commit deb4cc6

Please sign in to comment.