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/hpolytope.h b/include/convex_bodies/hpolytope.h
index d0e894faa..b01b6307c 100644
--- a/include/convex_bodies/hpolytope.h
+++ b/include/convex_bodies/hpolytope.h
@@ -6,6 +6,7 @@
 //Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018-19 programs.
 //Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program.
 //Contributed and/or modified by Alexandros Manochis, as part of Google Summer of Code 2020 program.
+//Contributed and/or modified by Luca Perju, as part of Google Summer of Code 2024 program.
 
 // Licensed under GNU LGPL.3, see LICENCE file
 
@@ -56,6 +57,7 @@ class HPolytope {
     MT                   A; //matrix A
     VT                   b; // vector b, s.t.: Ax<=b
     std::pair<Point, NT> _inner_ball;
+    bool                 normalized = false; // true if the polytope is normalized
 
 public:
     //TODO: the default implementation of the Big3 should be ok. Recheck.
@@ -68,7 +70,7 @@ class HPolytope {
 
     // Copy constructor
     HPolytope(HPolytope<Point> const& p) :
-            _d{p._d}, A{p.A}, b{p.b}, _inner_ball{p._inner_ball}
+            _d{p._d}, A{p.A}, b{p.b}, _inner_ball{p._inner_ball}, normalized{p.normalized}
     {
     }
 
@@ -172,11 +174,17 @@ class HPolytope {
         return b;
     }
 
+    bool is_normalized ()
+    {
+        return normalized;
+    }
+
 
     // change the matrix A
     void set_mat(MT const& A2)
     {
         A = A2;
+        normalized = false;
     }
 
 
@@ -541,7 +549,6 @@ class HPolytope {
 
         NT lamda = 0;
         VT sum_nom;
-        unsigned int j;
         int m = num_of_hyperplanes(), facet;
 
         Ar.noalias() += lambda_prev*Av;
@@ -822,6 +829,7 @@ class HPolytope {
     void linear_transformIt(MT const& T)
     {
         A = A * T;
+        normalized = false;
     }
 
 
@@ -866,6 +874,7 @@ class HPolytope {
             A.row(i) = A.row(i) / row_norm;
             b(i) = b(i) / row_norm;
         }
+        normalized = true;
     }
 
     void compute_reflection(Point& v, Point const&, int const& facet) const
@@ -909,6 +918,18 @@ class HPolytope {
             v += a;
     }
 
+    template <typename update_parameters>
+    NT compute_reflection(Point &v, const Point &, MT const &AE, VT const &AEA, NT const &vEv, update_parameters const &params) const {
+
+            Point a((-2.0 * params.inner_vi_ak) * A.row(params.facet_prev));
+            VT x = v.getCoefficients();
+            NT new_vEv = vEv - (4.0 * params.inner_vi_ak) * (AE.row(params.facet_prev).dot(x) - params.inner_vi_ak * AEA(params.facet_prev));
+            v += a;
+            NT coeff = std::sqrt(vEv / new_vEv);
+            v = v * coeff;
+            return coeff;
+    }
+
     void update_position_internal(NT&){}
 
     template <class bfunc, class NonLinearOracle>
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 <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,6 +53,23 @@ 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
@@ -45,7 +77,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 +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
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<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;
     }
 };
diff --git a/include/random_walks/gaussian_accelerated_billiard_walk.hpp b/include/random_walks/gaussian_accelerated_billiard_walk.hpp
index fbf43f443..d50cb7e5f 100644
--- a/include/random_walks/gaussian_accelerated_billiard_walk.hpp
+++ b/include/random_walks/gaussian_accelerated_billiard_walk.hpp
@@ -3,14 +3,18 @@
 // Copyright (c) 2012-2020 Vissarion Fisikopoulos
 // Copyright (c) 2018-2020 Apostolos Chalkis
 // Copyright (c) 2021 Vaibhav Thakkar
+// Copyright (c) 2024 Luca Perju
 
 // Contributed and/or modified by Vaibhav Thakkar, as part of Google Summer of Code 2021 program.
+// Contributed and/or modified by Luca Perju, as part of Google Summer of Code 2024 program.
 
 // Licensed under GNU LGPL.3, see LICENCE file
 
 #ifndef RANDOM_WALKS_GAUSSIAN_ACCELERATED_BILLIARD_WALK_HPP
 #define RANDOM_WALKS_GAUSSIAN_ACCELERATED_BILLIARD_WALK_HPP
 
+#include <Eigen/Eigen>
+#include <cmath>
 #include "convex_bodies/orderpolytope.h"
 #include "convex_bodies/ellipsoid.h"
 #include "convex_bodies/ballintersectconvex.h"
@@ -20,16 +24,14 @@
 
 #include "random_walks/compute_diameter.hpp"
 
-// Billiard walk which accelarates each step for uniform distribution and also takes into account
-// the shape of the polytope for generating directions.
 
-struct GaussianAcceleratedBilliardWalk
+struct GABW
 {
-    GaussianAcceleratedBilliardWalk(double L)
+    GABW(double L)
             :   param(L, true)
     {}
 
-    GaussianAcceleratedBilliardWalk()
+    GABW()
             :   param(0, false)
     {}
 
@@ -57,56 +59,67 @@ struct GaussianAcceleratedBilliardWalk
 
 
     template
-            <
-                    typename Polytope,
-                    typename RandomNumberGenerator
-            >
+    <
+            typename Polytope,
+            typename RandomNumberGenerator
+    >
     struct Walk
     {
         typedef typename Polytope::PointType Point;
-        // typedef typename Polytope::MT MT;
+        typedef typename Polytope::MT MT;
+        typedef typename Polytope::VT VT;
         typedef typename Point::FT NT;
 
-        template <typename GenericPolytope, typename Ellipsoid>
+        template <typename GenericPolytope>
         Walk(GenericPolytope& P,
              Point const& p,
-             Ellipsoid const& E,   // ellipsoid representing the Gaussian distribution
+             MT const& E,   // covariance matrix representing the Gaussian distribution
              RandomNumberGenerator &rng)
         {
+            if(!P.is_normalized()) {
+                P.normalize();
+            }
             _update_parameters = update_parameters();
-            _Len = compute_diameter<GenericPolytope>
-                ::template compute<NT>(P);
-
-            // Removed as will be used for sparse matrices only
-            // _AA.noalias() = P.get_mat() * P.get_mat().transpose();
-            initialize(P, p, E, rng);
+            _L = compute_diameter<GenericPolytope>::template compute<NT>(P);
+            Eigen::LLT<MT> lltOfE(E.llt().solve(MT::Identity(E.cols(), E.cols()))); // compute the Cholesky decomposition of inv(E)
+            if (lltOfE.info() != Eigen::Success) {
+                throw std::runtime_error("Cholesky decomposition failed!");
+            }
+            _L_cov = lltOfE.matrixL();
+            _E = E;
+            _AA.noalias() = P.get_mat() * P.get_mat().transpose();
+            _rho = 1000 * P.dimension(); // upper bound for the number of reflections (experimental)
+            initialize(P, p, rng);
         }
 
-        template <typename GenericPolytope, typename Ellipsoid>
+        template <typename GenericPolytope>
         Walk(GenericPolytope& P,
              Point const& p,
-             Ellipsoid const& E,   // ellipsoid representing the Gaussian distribution
+             MT const& E,   // covariance matrix representing the Gaussian distribution
              RandomNumberGenerator &rng,
              parameters const& params)
         {
+            if(!P.is_normalized()) {
+                P.normalize();
+            }
             _update_parameters = update_parameters();
-            _Len = params.set_L ? params.m_L
+            _L = params.set_L ? params.m_L
                               : compute_diameter<GenericPolytope>
                                 ::template compute<NT>(P);
-
-            // Removed as will be used for sparse matrices only
-            // _AA.noalias() = P.get_mat() * P.get_mat().transpose();
-            initialize(P, p, E, rng);
+            Eigen::LLT<MT> lltOfE(E.llt().solve(MT::Identity(E.cols(), E.cols()))); // compute the Cholesky decomposition of inv(E)
+            if (lltOfE.info() != Eigen::Success) {
+                throw std::runtime_error("Cholesky decomposition failed!");
+            }
+            _L_cov = lltOfE.matrixL();
+            _E = E;
+            _AA.noalias() = P.get_mat() * P.get_mat().transpose();
+            _rho = 1000 * P.dimension(); // upper bound for the number of reflections (experimental)
+            initialize(P, p, rng);
         }
 
-        template
-                <
-                        typename GenericPolytope,
-                        typename Ellipsoid
-                >
+        template <typename GenericPolytope>
         inline void apply(GenericPolytope& P,
                           Point &p,       // a point to return the result
-                          Ellipsoid const& E,   // ellipsoid representing the Gaussian distribution
                           unsigned int const& walk_length,
                           RandomNumberGenerator &rng)
         {
@@ -114,25 +127,50 @@ struct GaussianAcceleratedBilliardWalk
             NT T;
             const NT dl = 0.995;
             int it;
+            NT coef = 1.0;
+            NT vEv;
 
             for (auto j=0u; j<walk_length; ++j)
             {
-                T = -std::log(rng.sample_urdist()) * _Len;
-                _v = GetGaussianDirection<Point>::apply(n, E, rng);
-                NT norm_v = _v.length();
-                _v /= norm_v;
+                T = -std::log(rng.sample_urdist()) * _L;
+
+                _v = GetDirection<Point>::apply(n, rng, false);
+                _v = Point(_L_cov.template triangularView<Eigen::Lower>() * _v.getCoefficients());
+                coef = 1.0;
+
+                vEv = (_v.getCoefficients().transpose() * _E.template selfadjointView<Eigen::Upper>()).dot(_v.getCoefficients());
 
                 Point p0 = _p;
-                Point v0 = _v;
-                typename Point::Coeff lambdas0 = _lambdas;
-                typename Point::Coeff Av0 = _Av;
-                NT lambda_prev0 = _lambda_prev;
 
                 it = 0;
-                while (it < 100*n)
+                std::pair<NT, int> pbpair;
+                if(!was_reset) {
+                    pbpair = P.line_positive_intersect(_p, _v, _lambdas, _Av, _lambda_prev, _update_parameters);
+                } else {
+                    pbpair = P.line_first_positive_intersect(_p, _v, _lambdas, _Av, _update_parameters);
+                    was_reset = false;
+                }                                                    
+
+                if (T <= pbpair.first) {
+                    _p += (T * _v);
+                    _lambda_prev = T;
+                    continue;
+                }
+
+                _lambda_prev = dl * pbpair.first;
+                _p += (_lambda_prev * _v);
+                T -= _lambda_prev;
+                coef = P.compute_reflection(_v, _p, _AE, _AEA, vEv, _update_parameters);
+                it++;
+
+                while (it < _rho)
                 {
                     std::pair<NT, int> pbpair
-                            = P.line_positive_intersect(_p, _v, _lambdas, _Av, _lambda_prev, _update_parameters);
+                        = P.line_positive_intersect(_p, _v, _lambdas, _Av, _lambda_prev, _AA, _update_parameters);
+                    _Av *= coef;
+                    _update_parameters.inner_vi_ak *= coef;
+                    pbpair.first /= coef;
+
                     if (T <= pbpair.first) {
                         _p += (T * _v);
                         _lambda_prev = T;
@@ -141,63 +179,42 @@ struct GaussianAcceleratedBilliardWalk
                     _lambda_prev = dl * pbpair.first;
                     _p += (_lambda_prev * _v);
                     T -= _lambda_prev;
-                    P.compute_reflection(_v, _p, _update_parameters);
+                    coef = P.compute_reflection(_v, _p, _AE, _AEA, vEv, _update_parameters);
                     it++;
                 }
-                if (it == 100*n)
-                {
+                if (it == _rho){
                     _p = p0;
-                    _lambdas = lambdas0;
-                    _Av = Av0;
-                    _lambda_prev = lambda_prev0;
-                }
-                else
-                {
-                    // metropolis filter
-                    NT u_logprob = log(rng.sample_urdist());
-                    NT accept_logprob = 0.5 * (norm_v * norm_v) * (E.mat_mult(v0) - E.mat_mult(_v));
-                    // std::cout << "diff: " << (accept_logprob - u_logprob) << std::endl;
-                    if(u_logprob > accept_logprob) {
-                        // reject
-                        _p = p0;
-                        _lambdas = lambdas0;
-                        _Av = Av0;
-                        _lambda_prev = lambda_prev0;
-                    }
+                    was_reset = true;
                 }
             }
-            p = _p;
-        }
 
-        inline void update_delta(NT L)
-        {
-            _Len = L;
+            p = _p;
         }
 
     private :
 
-        template
-                <
-                        typename GenericPolytope,
-                        typename Ellipsoid
-                >
+        template <typename GenericPolytope>
         inline void initialize(GenericPolytope& P,
                                Point const& p,  // a point to start
-                               Ellipsoid const& E,   // ellipsoid representing the Gaussian distribution
                                RandomNumberGenerator &rng)
         {
             unsigned int n = P.dimension();
             const NT dl = 0.995;
+            was_reset = false;
             _lambdas.setZero(P.num_of_hyperplanes());
             _Av.setZero(P.num_of_hyperplanes());
             _p = p;
-            _v = GetGaussianDirection<Point>::apply(n, E, rng);
-            NT norm_v = _v.length();
-            _v /= norm_v;
+            _AE.noalias() = P.get_mat() * _E;
+            _AEA = _AE.cwiseProduct(P.get_mat()).rowwise().sum();
+
+            _v = GetDirection<Point>::apply(n, rng, false);
+            _v = Point(_L_cov.template triangularView<Eigen::Lower>() * _v.getCoefficients());
 
-            NT T = -std::log(rng.sample_urdist()) * _Len;
+            NT T = -std::log(rng.sample_urdist()) * _L;
             Point p0 = _p;
             int it = 0;
+            NT coef = 1.0;
+            NT vEv = (_v.getCoefficients().transpose() * _E.template selfadjointView<Eigen::Upper>()).dot(_v.getCoefficients());
 
             std::pair<NT, int> pbpair
                     = P.line_first_positive_intersect(_p, _v, _lambdas, _Av, _update_parameters);
@@ -209,17 +226,21 @@ struct GaussianAcceleratedBilliardWalk
             _lambda_prev = dl * pbpair.first;
             _p += (_lambda_prev * _v);
             T -= _lambda_prev;
-            P.compute_reflection(_v, _p, _update_parameters);
+            coef = P.compute_reflection(_v, _p, _AE, _AEA, vEv, _update_parameters);
 
-            while (it <= 100*n)
+            while (it <= _rho)
             {
                 std::pair<NT, int> pbpair
-                        = P.line_positive_intersect(_p, _v, _lambdas, _Av, _lambda_prev, _update_parameters);
+                        = P.line_positive_intersect(_p, _v, _lambdas, _Av, _lambda_prev, _AA, _update_parameters);
+                _Av *= coef;
+                _update_parameters.inner_vi_ak *= coef;
+                pbpair.first /= coef;
+
                 if (T <= pbpair.first) {
                     _p += (T * _v);
                     _lambda_prev = T;
                     break;
-                } else if (it == 100*n) {
+                } else if (it == _rho) {
                     _lambda_prev = rng.sample_urdist() * pbpair.first;
                     _p += (_lambda_prev * _v);
                     break;
@@ -227,19 +248,25 @@ struct GaussianAcceleratedBilliardWalk
                 _lambda_prev = dl * pbpair.first;
                 _p += (_lambda_prev * _v);
                 T -= _lambda_prev;
-                P.compute_reflection(_v, _p, _update_parameters);
+                coef = P.compute_reflection(_v, _p, _AE, _AEA, vEv, _update_parameters);
                 it++;
             }
         }
 
-        NT _Len;
+        NT _L;
         Point _p;
         Point _v;
         NT _lambda_prev;
-        // MT _AA;  // Removed as will be used for sparse matrices only
+        MT _AA;
+        MT _L_cov;   // LL' = inv(E)
+        MT _AE;
+        MT _E;
+        VT _AEA;
+        unsigned int _rho;
         update_parameters _update_parameters;
         typename Point::Coeff _lambdas;
         typename Point::Coeff _Av;
+        bool was_reset;
     };
 
 };
diff --git a/include/random_walks/uniform_accelerated_billiard_walk.hpp b/include/random_walks/uniform_accelerated_billiard_walk.hpp
index 980c8e5da..fa9c6fb8f 100644
--- a/include/random_walks/uniform_accelerated_billiard_walk.hpp
+++ b/include/random_walks/uniform_accelerated_billiard_walk.hpp
@@ -4,6 +4,7 @@
 // Copyright (c) 2018-2020 Apostolos Chalkis
 
 // Contributed and/or modified by Alexandros Manochis, as part of Google Summer of Code 2020 program.
+// Contributed and/or modified by Luca Perju, as part of Google Summer of Code 2024 program.
 
 // Licensed under GNU LGPL.3, see LICENCE file
 
@@ -49,10 +50,10 @@ struct AcceleratedBilliardWalk
 
 
     template
-            <
-                    typename Polytope,
-                    typename RandomNumberGenerator
-            >
+    <
+            typename Polytope,
+            typename RandomNumberGenerator
+    >
     struct Walk
     {
         typedef typename Polytope::PointType Point;
@@ -62,6 +63,9 @@ struct AcceleratedBilliardWalk
         template <typename GenericPolytope>
         Walk(GenericPolytope &P, Point const& p, RandomNumberGenerator &rng)
         {
+            if(!P.is_normalized()) {
+                P.normalize();
+            }
             _update_parameters = update_parameters();
             _L = compute_diameter<GenericPolytope>
                 ::template compute<NT>(P);
@@ -74,6 +78,9 @@ struct AcceleratedBilliardWalk
         Walk(GenericPolytope &P, Point const& p, RandomNumberGenerator &rng,
              parameters const& params)
         {
+            if(!P.is_normalized()) {
+                P.normalize();
+            }
             _update_parameters = update_parameters();
             _L = params.set_L ? params.m_L
                               : compute_diameter<GenericPolytope>
@@ -104,8 +111,14 @@ struct AcceleratedBilliardWalk
                 Point p0 = _p;
 
                 it = 0;
-                std::pair<NT, int> pbpair = P.line_positive_intersect(_p, _v, _lambdas, _Av,
-                                                                      _lambda_prev, _update_parameters);
+                std::pair<NT, int> pbpair;
+                if(!was_reset) {
+                    pbpair = P.line_positive_intersect(_p, _v, _lambdas, _Av, _lambda_prev, _update_parameters);
+                } else {
+                    pbpair = P.line_first_positive_intersect(_p, _v, _lambdas, _Av, _update_parameters);
+                    was_reset = false;
+                }
+                
                 if (T <= pbpair.first) {
                     _p += (T * _v);
                     _lambda_prev = T;
@@ -134,7 +147,10 @@ struct AcceleratedBilliardWalk
                     P.compute_reflection(_v, _p, _update_parameters);
                     it++;
                 }
-                if (it == _rho) _p = p0;
+                if (it == _rho) {
+                    _p = p0;
+                    was_reset = true;
+                }
             }
             p = _p;
         }
@@ -185,8 +201,7 @@ struct AcceleratedBilliardWalk
 
                 apply(P, p, walk_length, rng);
                 max_dist = get_max_distance(pointset, p, rad);
-                if (max_dist > Lmax) 
-                {
+                if (max_dist > Lmax) {
                     Lmax = max_dist;
                 }
                 if (2.0*rad > Lmax) {
@@ -196,8 +211,7 @@ struct AcceleratedBilliardWalk
             }
 
             if (Lmax > _L) {
-                if (P.dimension() <= 2500) 
-                {
+                if (P.dimension() <= 2500) {
                     update_delta(Lmax);
                 }
                 else{
@@ -231,6 +245,7 @@ struct AcceleratedBilliardWalk
         {
             unsigned int n = P.dimension();
             const NT dl = 0.995;
+            was_reset = false;
             _lambdas.setZero(P.num_of_hyperplanes());
             _Av.setZero(P.num_of_hyperplanes());
             _p = p;
@@ -301,6 +316,7 @@ struct AcceleratedBilliardWalk
         update_parameters _update_parameters;
         typename Point::Coeff _lambdas;
         typename Point::Coeff _Av;
+        bool was_reset;
     };
 
 };
diff --git a/include/random_walks/uniform_accelerated_billiard_walk_parallel.hpp b/include/random_walks/uniform_accelerated_billiard_walk_parallel.hpp
index 1c2630f9c..7ec437b23 100644
--- a/include/random_walks/uniform_accelerated_billiard_walk_parallel.hpp
+++ b/include/random_walks/uniform_accelerated_billiard_walk_parallel.hpp
@@ -53,10 +53,10 @@ struct AcceleratedBilliardWalkParallel
 
 
     template
-            <
-                    typename Polytope,
-                    typename RandomNumberGenerator
-            >
+    <
+            typename Polytope,
+            typename RandomNumberGenerator
+    >
     struct Walk
     {
         typedef typename Polytope::PointType Point;
@@ -66,6 +66,9 @@ struct AcceleratedBilliardWalkParallel
         template <typename GenericPolytope>
         Walk(GenericPolytope &P)
         {
+            if(!P.is_normalized()) {
+                P.normalize();
+            }
             _L = compute_diameter<GenericPolytope>
                 ::template compute<NT>(P);
             _AA.noalias() = P.get_mat() * P.get_mat().transpose();
@@ -76,6 +79,9 @@ struct AcceleratedBilliardWalkParallel
         template <typename GenericPolytope>
         Walk(GenericPolytope &P, NT const& L)
         {
+            if(!P.is_normalized()) {
+                P.normalize();
+            }
             _L = L > NT(0) ? L
                               : compute_diameter<GenericPolytope>
                                 ::template compute<NT>(P);
@@ -190,8 +196,7 @@ struct AcceleratedBilliardWalkParallel
 
                 apply(P, params, walk_length, rng);
                 max_dist = get_max_distance(pointset, params.p, rad);
-                if (max_dist > Lmax) 
-                {
+                if (max_dist > Lmax) {
                     Lmax = max_dist;
                 }
                 if (2.0*rad > Lmax) {
@@ -200,10 +205,8 @@ struct AcceleratedBilliardWalkParallel
                 pointset.push_back(params.p);
             }
 
-            if (Lmax > _L) 
-            {
-                if (P.dimension() <= 2500) 
-                {
+            if (Lmax > _L) {
+                if (P.dimension() <= 2500) {
                     update_delta(Lmax);
                 }
                 else{
diff --git a/include/sampling/random_point_generators.hpp b/include/sampling/random_point_generators.hpp
index e70e14617..e38f50cac 100644
--- a/include/sampling/random_point_generators.hpp
+++ b/include/sampling/random_point_generators.hpp
@@ -121,7 +121,7 @@ struct MultivariateGaussianRandomPointGenerator
         Walk walk(P, p, E, rng);
         for (unsigned int i=0; i<rnum; ++i)
         {
-            walk.template apply(P, p, E, walk_length, rng);
+            walk.template apply(P, p, walk_length, rng);
             policy.apply(randPoints, p);
         }
     }
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 615791148..61e499a9e 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -209,6 +209,7 @@ add_test(NAME test_bcdhr COMMAND sampling_test -tc=bcdhr)
 add_test(NAME test_grdhr COMMAND sampling_test -tc=grdhr)
 add_test(NAME test_gbaw COMMAND sampling_test -tc=gbaw)
 add_test(NAME test_ghmc COMMAND sampling_test -tc=ghmc)
+add_test(NAME test_gabw COMMAND sampling_test -tc=gabw)
 
 add_executable (mmcs_test mmcs_test.cpp $<TARGET_OBJECTS:test_main>)
 add_test(NAME test_mmcs COMMAND mmcs_test -tc=mmcs)
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);
+    
 }
 
 
diff --git a/test/sampling_test.cpp b/test/sampling_test.cpp
index 8e82e5818..1929c25fd 100644
--- a/test/sampling_test.cpp
+++ b/test/sampling_test.cpp
@@ -28,6 +28,8 @@
 
 #include "diagnostics/univariate_psrf.hpp"
 
+#include "preprocess/inscribed_ellipsoid_rounding.hpp"
+
 
 template
 <
@@ -302,6 +304,53 @@ void call_test_ghmc(){
     CHECK(score.maxCoeff() < 2.2);
 }
 
+template <typename NT, typename WalkType = GABW>
+void call_test_gabw(){
+    typedef Cartesian<NT>    Kernel;
+    typedef typename Kernel::Point    Point;
+    typedef HPolytope<Point> Hpolytope;
+    typedef Eigen::Matrix<NT,Eigen::Dynamic,Eigen::Dynamic> MT;
+    typedef Eigen::Matrix<NT,Eigen::Dynamic,1> VT;
+    Hpolytope P;
+    unsigned int d = 10;
+    typedef BoostRandomNumberGenerator<boost::mt19937, NT, 3> RNGType;
+
+    unsigned int walkL = 10, numpoints = 10000, nburns = 0;
+    RNGType rng(d);
+    Point StartingPoint(d);
+    std::list<Point> randPoints;
+
+    std::cout << "--- Testing Gaussian Accelerated Billiard Walk for Skinny-H-cube10" << std::endl;
+    P = generate_skinny_cube<Hpolytope>(10);
+
+
+    Point p = P.ComputeInnerBall().first;
+    typedef typename GABW::template Walk
+            <
+                    Hpolytope,
+                    RNGType
+            > walk;
+    typedef MultivariateGaussianRandomPointGenerator <walk> RandomPointGenerator;
+    PushBackWalkPolicy push_back_policy;
+
+    std::tuple<MT, VT, NT> ellipsoid = compute_inscribed_ellipsoid<MT, EllipsoidType::MAX_ELLIPSOID>
+    (P.get_mat(), P.get_vec(), p.getCoefficients(), 500, std::pow(10, -6.0), std::pow(10, -4.0));
+    const MT E = get<0>(ellipsoid);
+
+    RandomPointGenerator::apply(P, p, E, numpoints, 1, randPoints,
+                                push_back_policy, rng);
+
+    MT samples(d, numpoints);
+    unsigned int jj = 0;
+    for (typename std::list<Point>::iterator rpit = randPoints.begin(); rpit!=randPoints.end(); rpit++, jj++)
+        samples.col(jj) = (*rpit).getCoefficients();
+
+    VT score = univariate_psrf<NT, VT>(samples);
+    std::cout << "psrf = " << score.maxCoeff() << std::endl;
+
+    CHECK(score.maxCoeff() < 1.1);
+}
+
 TEST_CASE("dikin") {
     call_test_dikin<double>();
 }
@@ -333,3 +382,7 @@ TEST_CASE("gbaw") {
 TEST_CASE("ghmc") {
     call_test_ghmc<double>();
 }
+
+TEST_CASE("gabw") {
+    call_test_gabw<double>();
+}