From b2806e47738614abd7d51ccf8ded797bdf0fcc7e Mon Sep 17 00:00:00 2001
From: lucaperju <lucaperju@gmail.com>
Date: Fri, 28 Jun 2024 17:20:00 +0200
Subject: [PATCH 1/6] Order Polytopes generation

---
 examples/crhmc_sampling/.gitignore            |  1 +
 examples/crhmc_sampling/simple_crhmc.cpp      | 96 +++++++++++++++----
 include/convex_bodies/orderpolytope.h         |  6 ++
 include/generators/order_polytope_generator.h | 79 ++++++++++++++-
 include/misc/poset.h                          | 65 +++++++------
 include/preprocess/svd_rounding.hpp           |  4 +-
 6 files changed, 198 insertions(+), 53 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/examples/crhmc_sampling/simple_crhmc.cpp b/examples/crhmc_sampling/simple_crhmc.cpp
index 8f07065b0..733d2031e 100644
--- a/examples/crhmc_sampling/simple_crhmc.cpp
+++ b/examples/crhmc_sampling/simple_crhmc.cpp
@@ -15,6 +15,8 @@
 // Monte Carlo in a Constrained Space"
 #include "Eigen/Eigen"
 #include "cartesian_geom/cartesian_kernel.h"
+#include "generators/h_polytopes_generator.h"
+#include "generators/order_polytope_generator.h"
 #include "volume/sampling_policies.hpp"
 #include "ode_solvers/ode_solvers.hpp"
 #include "preprocess/crhmc/crhmc_input.h"
@@ -22,15 +24,42 @@
 #include "sampling/random_point_generators.hpp"
 #include "sampling/sampling.hpp"
 #include "misc/misc.h"
+#include "misc/poset.h"
 #include "random.hpp"
 #include <vector>
 #include "random_walks/random_walks.hpp"
 #include "generators/known_polytope_generators.h"
 #include "helper_functions.hpp"
+#include "volume/rotating.hpp"
+#include "convex_bodies/orderpolytope.h"
+
+
+template 
+<
+typename Polytope,
+typename Point,
+typename RandomNumberGenerator
+>
+void sample_cdhr (Polytope &P, RandomNumberGenerator &rng, std::list<Point> &randPoints, unsigned int const&N)
+{
+        Point p = P.ComputeInnerBall().first;
+        typedef typename CDHRWalk::template Walk
+                <
+                        Polytope,
+                        RandomNumberGenerator
+                > walk;
+
+        typedef RandomPointGenerator <walk> RandomPointGenerator;
+        PushBackWalkPolicy push_back_policy;
+
+        RandomPointGenerator::apply(P, p, N, 1, randPoints,
+                                    push_back_policy, rng);
+}
+
 
 template <int simdLen>
 void sample_hpoly(int n_samples = 80000,
-              int n_burns = 20000, int dim = 2) {
+              int n_burns = 20000, int dim = 20, int m = 150, bool order_poly = false, bool crhmc_walk = false) {
   using NT = double;
   using Kernel = Cartesian<NT>;
   using Point = typename Kernel::Point;
@@ -40,18 +69,46 @@ void sample_hpoly(int n_samples = 80000,
   using PolytopeType = HPolytope<Point>;
   using VT = Eigen::Matrix<NT, Eigen::Dynamic, 1>;
   using MT = PolytopeType::MT;
+  typedef boost::mt19937 PolyRNGType;
   using RNG = BoostRandomNumberGenerator<boost::mt19937, NT>;
-  std::string problem_name("simplex");
-  std::cerr << "CRHMC on " << problem_name << "\n";
-  RNG rng(1);
-  PolytopeType HP=generate_simplex<PolytopeType>(dim,false);
-  int dimension = HP.dimension();
+
+
+  RNG rng(dim);
+  PolytopeType HP;
+  if(order_poly) {
+    HP = random_orderpoly<PolytopeType, NT>(dim, m);
+    std::cout << "Sampling from Order Polytope" << std::endl;
+  } else {
+    HP = skinny_random_hpoly<PolytopeType, NT, PolyRNGType>(dim, m, false, NT(4000));
+    std::cout << "Sampling from Random skinny Polytope" << std::endl;
+    // HP = generate_skinny_cube<PolytopeType>(20);
+    // rotating<MT>(HP);
+  }
+
+  // HP.print();
   Func * f = new Func;
   Grad * g = new Grad;
   std::list<Point> PointList;
-  execute_crhmc< PolytopeType, RNG, std::list<Point>, Grad, Func, Hess, CRHMCWalk, simdLen>(
+
+  
+  std::chrono::time_point<std::chrono::high_resolution_clock> start, stop;
+  start = std::chrono::high_resolution_clock::now();
+
+  if(crhmc_walk) {
+    std::cout << "Using CRHMC walk" << std::endl;
+    execute_crhmc< PolytopeType, RNG, std::list<Point>, Grad, Func, Hess, CRHMCWalk, simdLen>(
       HP, rng, PointList, 1, n_samples, n_burns, g, f);
-  MT samples = MT(dimension, PointList.size());
+  } else {
+    std::cout << "Using CDHR walk" << std::endl;
+    sample_cdhr(HP, rng, PointList, n_samples);
+  }
+
+  stop = std::chrono::high_resolution_clock::now();
+
+  std::chrono::duration<double> total_time = stop - start;
+  std::cout << "Done in " << total_time.count() << '\n';
+
+  MT samples = MT(dim, PointList.size());
   int i=0;
   for (std::list<Point>::iterator it = PointList.begin(); it != PointList.end(); ++it){
     samples.col(i) = (*it).getCoefficients();
@@ -59,36 +116,35 @@ void sample_hpoly(int n_samples = 80000,
   }
   std::cerr<<"max_psrf: "<< max_interval_psrf<NT,VT,MT>(samples)<<"\n";
   std::ofstream samples_stream;
-  samples_stream.open("CRHMC_SIMD_" + std::to_string(simdLen) + "_" +
-                      problem_name + "_samples.txt");
+  samples_stream.open("CRHMC_SIMD_" + std::to_string(simdLen) + "_simplex" + "_samples.txt");
   samples_stream << samples.transpose() << std::endl;
   delete f;
   delete g;
 }
 
+
 template<int simdLen>
 void run_main(int n_samples = 80000,
               int n_burns = 20000,
-              int dimension = 2){
-  std::cerr<<"Sampling HPolytope\n";
-  sample_hpoly<simdLen>(n_samples, n_burns, dimension);
+              int dimension = 20, int m = 150, bool order_poly = false, bool crhmc_walk = false){
+  sample_hpoly<simdLen>(n_samples, n_burns, dimension, m, order_poly, crhmc_walk);
 }
 int main(int argc, char *argv[]) {
-  if (argc != 5) {
+  if (argc != 8) {
     std::cerr << "Example Usage: ./simple_crhmc "
-                 "[simdLen] [n_samples] [n_burns] [dimension]\n";
-    std::cerr << "i.e.: ./simple_crhmc 4 1000 500 2\n";
+                 "[simdLen] [n_samples] [n_burns] [dimension] [facets] [if_order_poly] [if_crhmc_walk]\n";
+    std::cerr << "i.e.: ./simple_crhmc 4 1000 500 20 150 1 0\n";
     exit(1);
   }
   std::cerr << "To plot: python3 ../python_utilities/plot_samples.py <CRHMC_SIMD_4_simplex_samples.txt --save"<<"\n";
   if (atoi(argv[1]) == 1) {
-    run_main<1>(atoi(argv[2]), atoi(argv[3]), atoi(argv[4]));
+    run_main<1>(atoi(argv[2]), atoi(argv[3]), atoi(argv[4]), atoi(argv[5]), atoi(argv[6]), atoi(argv[7]));
   } else if (atoi(argv[1]) == 4) {
-    run_main<4>(atoi(argv[2]), atoi(argv[3]), atoi(argv[4]));
+    run_main<4>(atoi(argv[2]), atoi(argv[3]), atoi(argv[4]), atoi(argv[5]), atoi(argv[6]), atoi(argv[7]));
   } else if (atoi(argv[1]) == 8) {
-    run_main<8>(atoi(argv[2]), atoi(argv[3]), atoi(argv[4]));
+    run_main<8>(atoi(argv[2]), atoi(argv[3]), atoi(argv[4]), atoi(argv[5]), atoi(argv[6]), atoi(argv[7]));
   } else if (atoi(argv[1]) == 16) {
-    run_main<16>(atoi(argv[2]), atoi(argv[3]), atoi(argv[4]));
+    run_main<16>(atoi(argv[2]), atoi(argv[3]), atoi(argv[4]), atoi(argv[5]), atoi(argv[6]), atoi(argv[7]));
   }
   return 0;
 }
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..6d8449c89 100644
--- a/include/generators/order_polytope_generator.h
+++ b/include/generators/order_polytope_generator.h
@@ -9,11 +9,25 @@
 #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"
+
 
 // Instances taken from: https://github.com/ttalvitie/le-counting-practice
 static const std::unordered_map<std::string, std::string> instances =
@@ -38,6 +52,21 @@ 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 &poset) {
+    typedef typename Polytope::PointType Point;
+
+    OrderPolytope<Point> OP(poset);
+    if constexpr (std::is_same< Polytope, OrderPolytope<Point> >::value ) {
+        return OP;
+    } else {
+        Polytope HP(OP.dimension(), OP.get_full_mat(), OP.get_vec());
+        return HP;
+    }
+}
+
 // 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 +74,7 @@ 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 +85,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)) {
+        int 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
diff --git a/include/misc/poset.h b/include/misc/poset.h
index 69ba02ec7..c8ddb94be 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(std::vector<unsigned int> &res, const unsigned int &n, const RV &relations)
+    {
+        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(order, n, relations);
+        
+        if(order.size() < n) {
+            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(res, n, order_relations);
         return res;
     }
 };
diff --git a/include/preprocess/svd_rounding.hpp b/include/preprocess/svd_rounding.hpp
index e62c03d28..ac8793e75 100644
--- a/include/preprocess/svd_rounding.hpp
+++ b/include/preprocess/svd_rounding.hpp
@@ -46,6 +46,8 @@ void svd_on_sample(Polytope &P, Point &p, unsigned int const& num_rounding_steps
         RetMat.row(jj) = (*rpit).getCoefficients().transpose();
     }
 
+    //std::cout << RetMat << std::endl;
+
     for (int i = 0; i < P.dimension(); ++i) {
         Means(i) = RetMat.col(i).mean();
     }
@@ -101,7 +103,7 @@ std::tuple<MT, VT, NT> svd_rounding(Polytope &P,
 
     bool done = false, last_round_under_p, fail;
 
-    unsigned int tries=0, num_rounding_steps = 10 * n, rounding_samples = 0, round_it;
+    unsigned int tries=0, num_rounding_steps = 50 * n, rounding_samples = 0, round_it;
     NT max_s, s_cutof, p_cutof, num_its, prev_max_s = std::numeric_limits<NT>::max(),
        s_cutoff, p_cutoff;
     MT V(n,n), S(n,n);

From c7c265b7818affd00ec5f77d257ed5c09f20e77f Mon Sep 17 00:00:00 2001
From: lucaperju <lucaperju@gmail.com>
Date: Thu, 11 Jul 2024 14:46:14 +0200
Subject: [PATCH 2/6] minor changes for PR

---
 examples/crhmc_sampling/simple_crhmc.cpp | 98 +++++-------------------
 1 file changed, 21 insertions(+), 77 deletions(-)

diff --git a/examples/crhmc_sampling/simple_crhmc.cpp b/examples/crhmc_sampling/simple_crhmc.cpp
index 733d2031e..0a68dbc33 100644
--- a/examples/crhmc_sampling/simple_crhmc.cpp
+++ b/examples/crhmc_sampling/simple_crhmc.cpp
@@ -15,8 +15,6 @@
 // Monte Carlo in a Constrained Space"
 #include "Eigen/Eigen"
 #include "cartesian_geom/cartesian_kernel.h"
-#include "generators/h_polytopes_generator.h"
-#include "generators/order_polytope_generator.h"
 #include "volume/sampling_policies.hpp"
 #include "ode_solvers/ode_solvers.hpp"
 #include "preprocess/crhmc/crhmc_input.h"
@@ -24,42 +22,15 @@
 #include "sampling/random_point_generators.hpp"
 #include "sampling/sampling.hpp"
 #include "misc/misc.h"
-#include "misc/poset.h"
 #include "random.hpp"
 #include <vector>
 #include "random_walks/random_walks.hpp"
 #include "generators/known_polytope_generators.h"
 #include "helper_functions.hpp"
-#include "volume/rotating.hpp"
-#include "convex_bodies/orderpolytope.h"
-
-
-template 
-<
-typename Polytope,
-typename Point,
-typename RandomNumberGenerator
->
-void sample_cdhr (Polytope &P, RandomNumberGenerator &rng, std::list<Point> &randPoints, unsigned int const&N)
-{
-        Point p = P.ComputeInnerBall().first;
-        typedef typename CDHRWalk::template Walk
-                <
-                        Polytope,
-                        RandomNumberGenerator
-                > walk;
-
-        typedef RandomPointGenerator <walk> RandomPointGenerator;
-        PushBackWalkPolicy push_back_policy;
-
-        RandomPointGenerator::apply(P, p, N, 1, randPoints,
-                                    push_back_policy, rng);
-}
-
 
 template <int simdLen>
 void sample_hpoly(int n_samples = 80000,
-              int n_burns = 20000, int dim = 20, int m = 150, bool order_poly = false, bool crhmc_walk = false) {
+              int n_burns = 20000, int dim = 2) {
   using NT = double;
   using Kernel = Cartesian<NT>;
   using Point = typename Kernel::Point;
@@ -69,46 +40,18 @@ void sample_hpoly(int n_samples = 80000,
   using PolytopeType = HPolytope<Point>;
   using VT = Eigen::Matrix<NT, Eigen::Dynamic, 1>;
   using MT = PolytopeType::MT;
-  typedef boost::mt19937 PolyRNGType;
   using RNG = BoostRandomNumberGenerator<boost::mt19937, NT>;
-
-
-  RNG rng(dim);
-  PolytopeType HP;
-  if(order_poly) {
-    HP = random_orderpoly<PolytopeType, NT>(dim, m);
-    std::cout << "Sampling from Order Polytope" << std::endl;
-  } else {
-    HP = skinny_random_hpoly<PolytopeType, NT, PolyRNGType>(dim, m, false, NT(4000));
-    std::cout << "Sampling from Random skinny Polytope" << std::endl;
-    // HP = generate_skinny_cube<PolytopeType>(20);
-    // rotating<MT>(HP);
-  }
-
-  // HP.print();
+  std::string problem_name("simplex");
+  std::cerr << "CRHMC on " << problem_name << "\n";
+  RNG rng(1);
+  PolytopeType HP=generate_simplex<PolytopeType>(dim,false);
+  int dimension = HP.dimension();
   Func * f = new Func;
   Grad * g = new Grad;
   std::list<Point> PointList;
-
-  
-  std::chrono::time_point<std::chrono::high_resolution_clock> start, stop;
-  start = std::chrono::high_resolution_clock::now();
-
-  if(crhmc_walk) {
-    std::cout << "Using CRHMC walk" << std::endl;
-    execute_crhmc< PolytopeType, RNG, std::list<Point>, Grad, Func, Hess, CRHMCWalk, simdLen>(
+  execute_crhmc< PolytopeType, RNG, std::list<Point>, Grad, Func, Hess, CRHMCWalk, simdLen>(
       HP, rng, PointList, 1, n_samples, n_burns, g, f);
-  } else {
-    std::cout << "Using CDHR walk" << std::endl;
-    sample_cdhr(HP, rng, PointList, n_samples);
-  }
-
-  stop = std::chrono::high_resolution_clock::now();
-
-  std::chrono::duration<double> total_time = stop - start;
-  std::cout << "Done in " << total_time.count() << '\n';
-
-  MT samples = MT(dim, PointList.size());
+  MT samples = MT(dimension, PointList.size());
   int i=0;
   for (std::list<Point>::iterator it = PointList.begin(); it != PointList.end(); ++it){
     samples.col(i) = (*it).getCoefficients();
@@ -116,35 +59,36 @@ void sample_hpoly(int n_samples = 80000,
   }
   std::cerr<<"max_psrf: "<< max_interval_psrf<NT,VT,MT>(samples)<<"\n";
   std::ofstream samples_stream;
-  samples_stream.open("CRHMC_SIMD_" + std::to_string(simdLen) + "_simplex" + "_samples.txt");
+  samples_stream.open("CRHMC_SIMD_" + std::to_string(simdLen) + "_" +
+                      problem_name + "_samples.txt");
   samples_stream << samples.transpose() << std::endl;
   delete f;
   delete g;
 }
 
-
 template<int simdLen>
 void run_main(int n_samples = 80000,
               int n_burns = 20000,
-              int dimension = 20, int m = 150, bool order_poly = false, bool crhmc_walk = false){
-  sample_hpoly<simdLen>(n_samples, n_burns, dimension, m, order_poly, crhmc_walk);
+              int dimension = 2){
+  std::cerr<<"Sampling HPolytope\n";
+  sample_hpoly<simdLen>(n_samples, n_burns, dimension);
 }
 int main(int argc, char *argv[]) {
-  if (argc != 8) {
+  if (argc != 5) {
     std::cerr << "Example Usage: ./simple_crhmc "
-                 "[simdLen] [n_samples] [n_burns] [dimension] [facets] [if_order_poly] [if_crhmc_walk]\n";
-    std::cerr << "i.e.: ./simple_crhmc 4 1000 500 20 150 1 0\n";
+                 "[simdLen] [n_samples] [n_burns] [dimension]\n";
+    std::cerr << "i.e.: ./simple_crhmc 4 1000 500 2\n";
     exit(1);
   }
   std::cerr << "To plot: python3 ../python_utilities/plot_samples.py <CRHMC_SIMD_4_simplex_samples.txt --save"<<"\n";
   if (atoi(argv[1]) == 1) {
-    run_main<1>(atoi(argv[2]), atoi(argv[3]), atoi(argv[4]), atoi(argv[5]), atoi(argv[6]), atoi(argv[7]));
+    run_main<1>(atoi(argv[2]), atoi(argv[3]), atoi(argv[4]));
   } else if (atoi(argv[1]) == 4) {
-    run_main<4>(atoi(argv[2]), atoi(argv[3]), atoi(argv[4]), atoi(argv[5]), atoi(argv[6]), atoi(argv[7]));
+    run_main<4>(atoi(argv[2]), atoi(argv[3]), atoi(argv[4]));
   } else if (atoi(argv[1]) == 8) {
-    run_main<8>(atoi(argv[2]), atoi(argv[3]), atoi(argv[4]), atoi(argv[5]), atoi(argv[6]), atoi(argv[7]));
+    run_main<8>(atoi(argv[2]), atoi(argv[3]), atoi(argv[4]));
   } else if (atoi(argv[1]) == 16) {
-    run_main<16>(atoi(argv[2]), atoi(argv[3]), atoi(argv[4]), atoi(argv[5]), atoi(argv[6]), atoi(argv[7]));
+    run_main<16>(atoi(argv[2]), atoi(argv[3]), atoi(argv[4]));
   }
   return 0;
-}
+}
\ No newline at end of file

From 11d04b6f47965f9ea79609f85b05c05c3d21b66c Mon Sep 17 00:00:00 2001
From: lucaperju <lucaperju@gmail.com>
Date: Fri, 12 Jul 2024 09:32:01 +0200
Subject: [PATCH 3/6] remove space and comment

---
 examples/crhmc_sampling/simple_crhmc.cpp | 2 +-
 include/preprocess/svd_rounding.hpp      | 4 +---
 2 files changed, 2 insertions(+), 4 deletions(-)

diff --git a/examples/crhmc_sampling/simple_crhmc.cpp b/examples/crhmc_sampling/simple_crhmc.cpp
index 0a68dbc33..8f07065b0 100644
--- a/examples/crhmc_sampling/simple_crhmc.cpp
+++ b/examples/crhmc_sampling/simple_crhmc.cpp
@@ -91,4 +91,4 @@ int main(int argc, char *argv[]) {
     run_main<16>(atoi(argv[2]), atoi(argv[3]), atoi(argv[4]));
   }
   return 0;
-}
\ No newline at end of file
+}
diff --git a/include/preprocess/svd_rounding.hpp b/include/preprocess/svd_rounding.hpp
index ac8793e75..e62c03d28 100644
--- a/include/preprocess/svd_rounding.hpp
+++ b/include/preprocess/svd_rounding.hpp
@@ -46,8 +46,6 @@ void svd_on_sample(Polytope &P, Point &p, unsigned int const& num_rounding_steps
         RetMat.row(jj) = (*rpit).getCoefficients().transpose();
     }
 
-    //std::cout << RetMat << std::endl;
-
     for (int i = 0; i < P.dimension(); ++i) {
         Means(i) = RetMat.col(i).mean();
     }
@@ -103,7 +101,7 @@ std::tuple<MT, VT, NT> svd_rounding(Polytope &P,
 
     bool done = false, last_round_under_p, fail;
 
-    unsigned int tries=0, num_rounding_steps = 50 * n, rounding_samples = 0, round_it;
+    unsigned int tries=0, num_rounding_steps = 10 * n, rounding_samples = 0, round_it;
     NT max_s, s_cutof, p_cutof, num_its, prev_max_s = std::numeric_limits<NT>::max(),
        s_cutoff, p_cutoff;
     MT V(n,n), S(n,n);

From 9b4d325c25b8171d61b09859b3725cee326f90d9 Mon Sep 17 00:00:00 2001
From: lucaperju <lucaperju@gmail.com>
Date: Fri, 12 Jul 2024 09:42:46 +0200
Subject: [PATCH 4/6] remove bug

---
 include/generators/order_polytope_generator.h | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/include/generators/order_polytope_generator.h b/include/generators/order_polytope_generator.h
index 6d8449c89..1cf020336 100644
--- a/include/generators/order_polytope_generator.h
+++ b/include/generators/order_polytope_generator.h
@@ -98,7 +98,7 @@ Polytope random_orderpoly(unsigned int dim, unsigned int m, int seed = std::nume
 
     int rng_seed = std::chrono::system_clock::now().time_since_epoch().count();
     if (!isnan(seed)) {
-        int rng_seed = seed;
+        rng_seed = seed;
     }
 
     typedef BoostRandomNumberGenerator<boost::mt19937, NT> RNG;

From 63300b25985b597a154ed7a06864d8fbe2b669c2 Mon Sep 17 00:00:00 2001
From: lucaperju <lucaperju@gmail.com>
Date: Fri, 12 Jul 2024 13:51:04 +0200
Subject: [PATCH 5/6] Unit test for Random Order Polytope, and minor changes

---
 include/generators/order_polytope_generator.h |  8 ++++++--
 include/misc/poset.h                          |  8 ++++----
 test/order_polytope.cpp                       | 14 ++++++++++++++
 3 files changed, 24 insertions(+), 6 deletions(-)

diff --git a/include/generators/order_polytope_generator.h b/include/generators/order_polytope_generator.h
index 1cf020336..e9199ba2c 100644
--- a/include/generators/order_polytope_generator.h
+++ b/include/generators/order_polytope_generator.h
@@ -27,6 +27,7 @@
 #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
@@ -55,15 +56,18 @@ 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 &poset) {
+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 {
+    } else if constexpr (std::is_same<Polytope, HPolytope<Point> >::value ){
         Polytope HP(OP.dimension(), OP.get_full_mat(), OP.get_vec());
         return HP;
+    } else {
+        // TODO: implement functionality for more polytope types
+        throw "Unable to generate an Order Polytope of requested type";
     }
 }
 
diff --git a/include/misc/poset.h b/include/misc/poset.h
index c8ddb94be..084bf66aa 100644
--- a/include/misc/poset.h
+++ b/include/misc/poset.h
@@ -25,7 +25,7 @@ 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(std::vector<unsigned int> &res, const unsigned int &n, const RV &relations)
+    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);
@@ -75,9 +75,9 @@ class Poset {
         }
 
         std::vector<unsigned int> order;
-        sorted_list(order, n, relations);
+        sorted_list(n, relations, order);
         
-        if(order.size() < n) {
+        if(order.size() < n) { // TODO: accept cycles in the poset
             throw "corresponding DAG is not acyclic";
         }
 
@@ -132,7 +132,7 @@ class Poset {
     std::vector<unsigned int> topologically_sorted_list() const
     {
         std::vector<unsigned int> res;
-        sorted_list(res, n, order_relations);
+        sorted_list(n, order_relations, res);
         return res;
     }
 };
diff --git a/test/order_polytope.cpp b/test/order_polytope.cpp
index e2ad4300d..9a83e13a4 100644
--- a/test/order_polytope.cpp
+++ b/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);
+    
 }
 
 

From 6e7c3f278edd5ab3ebfb98bd1c8e636d36746ccc Mon Sep 17 00:00:00 2001
From: lucaperju <lucaperju@gmail.com>
Date: Fri, 12 Jul 2024 18:05:39 +0200
Subject: [PATCH 6/6] remove comment

---
 include/generators/order_polytope_generator.h | 1 -
 1 file changed, 1 deletion(-)

diff --git a/include/generators/order_polytope_generator.h b/include/generators/order_polytope_generator.h
index e9199ba2c..3cd907bd9 100644
--- a/include/generators/order_polytope_generator.h
+++ b/include/generators/order_polytope_generator.h
@@ -66,7 +66,6 @@ Polytope get_orderpoly(Poset const &poset) {
         Polytope HP(OP.dimension(), OP.get_full_mat(), OP.get_vec());
         return HP;
     } else {
-        // TODO: implement functionality for more polytope types
         throw "Unable to generate an Order Polytope of requested type";
     }
 }