diff --git a/src/volesti/include/convex_bodies/ballintersectconvex.h b/src/volesti/include/convex_bodies/ballintersectconvex.h
index b61ae532..992aae1a 100644
--- a/src/volesti/include/convex_bodies/ballintersectconvex.h
+++ b/src/volesti/include/convex_bodies/ballintersectconvex.h
@@ -50,6 +50,12 @@ class BallIntersectPolytope {
         return P.get_vec();
     }
 
+    bool is_normalized () {
+        return true;
+    }
+
+    void normalize() {}
+
     int is_in(PointType const& p) const
     {
         if (B.is_in(p)==-1)
diff --git a/src/volesti/include/convex_bodies/correlation_matrices/correlation_spectrahedron.hpp b/src/volesti/include/convex_bodies/correlation_matrices/correlation_spectrahedron.hpp
index a615b33c..92c38dac 100755
--- a/src/volesti/include/convex_bodies/correlation_matrices/correlation_spectrahedron.hpp
+++ b/src/volesti/include/convex_bodies/correlation_matrices/correlation_spectrahedron.hpp
@@ -261,6 +261,12 @@ class CorrelationSpectrahedron : public Spectrahedron<Point>{
     MT get_mat() const {
         return MT::Identity(this->d, this->d);
     }
+
+    bool is_normalized () {
+        return true;
+    }
+
+    void normalize() {}
 };
 
 #endif //VOLESTI_CONVEX_BODIES_CORRELATION_MATRICES_VOLESTI_CORRELATION_SPECTRAHEDRON_HPP
\ No newline at end of file
diff --git a/src/volesti/include/convex_bodies/hpolytope.h b/src/volesti/include/convex_bodies/hpolytope.h
index d0e894fa..9cabcd9a 100644
--- a/src/volesti/include/convex_bodies/hpolytope.h
+++ b/src/volesti/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
 
@@ -40,22 +41,27 @@ bool is_inner_point_nan_inf(VT const& p)
 /// This class describes a polytope in H-representation or an H-polytope
 /// i.e. a polytope defined by a set of linear inequalities
 /// \tparam Point Point type
-template <typename Point>
+template 
+<
+    typename Point, 
+    typename MT_type = Eigen::Matrix<typename Point::FT, Eigen::Dynamic, Eigen::Dynamic>
+>
 class HPolytope {
 public:
     typedef Point                                             PointType;
     typedef typename Point::FT                                NT;
     typedef typename std::vector<NT>::iterator                viterator;
-    //using RowMatrixXd = Eigen::Matrix<NT, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
-    //typedef RowMatrixXd MT;
-    typedef Eigen::Matrix<NT, Eigen::Dynamic, Eigen::Dynamic> MT;
+    typedef MT_type                                           MT;
     typedef Eigen::Matrix<NT, Eigen::Dynamic, 1>              VT;
+    typedef Eigen::Matrix<NT, Eigen::Dynamic, Eigen::Dynamic> DenseMT;
 
 private:
     unsigned int         _d; //dimension
     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
+    bool                 has_ball = false;
 
 public:
     //TODO: the default implementation of the Big3 should be ok. Recheck.
@@ -66,9 +72,15 @@ class HPolytope {
     {
     }
 
+    template<typename T = DenseMT>
+    HPolytope(unsigned d_, DenseMT const& A_, VT const& b_, typename std::enable_if<!std::is_same<MT, T>::value, T>::type* = 0) :
+        _d{d_}, A{A_.sparseView()}, b{b_}
+    {
+    }
+
     // Copy constructor
-    HPolytope(HPolytope<Point> const& p) :
-            _d{p._d}, A{p.A}, b{p.b}, _inner_ball{p._inner_ball}
+    HPolytope(HPolytope<Point, MT> const& p) :
+            _d{p._d}, A{p.A}, b{p.b}, _inner_ball{p._inner_ball}, normalized{p.normalized}, has_ball{p.has_ball}
     {
     }
 
@@ -82,10 +94,10 @@ class HPolytope {
         for (unsigned int i = 1; i < Pin.size(); i++) {
             b(i - 1) = Pin[i][0];
             for (unsigned int j = 1; j < _d + 1; j++) {
-                A(i - 1, j - 1) = -Pin[i][j];
+                A.coeffRef(i - 1, j - 1) = -Pin[i][j];
             }
         }
-        _inner_ball.second = -1;
+        has_ball = false;
         //_inner_ball = ComputeChebychevBall<NT, Point>(A, b);
     }
 
@@ -98,41 +110,60 @@ class HPolytope {
     void set_InnerBall(std::pair<Point,NT> const& innerball) //const
     {
         _inner_ball = innerball;
+        has_ball = true;
     }
 
     void set_interior_point(Point const& r)
     {
         _inner_ball.first = r;
+        if(!is_in(r)) {
+            std::cerr << "point outside of polytope was set as interior point" << std::endl;
+            has_ball = false;
+            return;
+        }
+        if(is_normalized()) {
+            _inner_ball.second = (b - A * r.getCoefficients()).minCoeff();
+        } else {
+            _inner_ball.second = std::numeric_limits<NT>::max();
+            for(int i = 0; i < num_of_hyperplanes(); ++i) {
+                NT dist = (b(i) - A.row(i).dot(r.getCoefficients()) ) / A.row(i).norm();
+                if(dist < _inner_ball.second) {
+                    _inner_ball.second = dist;
+                }
+            }
+        }
+        has_ball = true;
     }
 
     //Compute Chebyshev ball of H-polytope P:= Ax<=b
-    //Use LpSolve library
+    //First try using max_inscribed_ball
+    //Use LpSolve library if it fails
     std::pair<Point, NT> ComputeInnerBall()
     {
         normalize();
-        #ifndef DISABLE_LPSOLVE
-            _inner_ball = ComputeChebychevBall<NT, Point>(A, b); // use lpsolve library
-        #else
-
-            if (_inner_ball.second <= NT(0)) {
-
-                NT const tol = 1e-08;
-                std::tuple<VT, NT, bool> inner_ball = max_inscribed_ball(A, b, 5000, tol);
-
-                // check if the solution is feasible
-                if (is_in(Point(std::get<0>(inner_ball))) == 0 || std::get<1>(inner_ball) < tol/2.0 ||
-                    std::isnan(std::get<1>(inner_ball)) || std::isinf(std::get<1>(inner_ball)) ||
-                    is_inner_point_nan_inf(std::get<0>(inner_ball)))
-                {
-                    _inner_ball.second = -1.0;
-                } else
-                {
-                    _inner_ball.first = Point(std::get<0>(inner_ball));
-                    _inner_ball.second = std::get<1>(inner_ball);
-                }
+        if (!has_ball) {
+            
+            has_ball = true;
+            NT const tol = 1e-08;
+            std::tuple<VT, NT, bool> inner_ball = max_inscribed_ball(A, b, 5000, tol);
+
+            // check if the solution is feasible
+            if (is_in(Point(std::get<0>(inner_ball))) == 0 || std::get<1>(inner_ball) < tol/2.0 ||
+                std::isnan(std::get<1>(inner_ball)) || std::isinf(std::get<1>(inner_ball)) ||
+                is_inner_point_nan_inf(std::get<0>(inner_ball))) {
+                
+                std::cerr << "Failed to compute max inscribed ball, trying to use lpsolve" << std::endl;
+                #ifndef DISABLE_LPSOLVE
+                    _inner_ball = ComputeChebychevBall<NT, Point>(A, b); // use lpsolve library
+                #else
+                    std::cerr << "lpsolve is disabled, unable to compute inner ball";
+                    has_ball = false;
+                #endif
+            } else {
+                _inner_ball.first = Point(std::get<0>(inner_ball));
+                _inner_ball.second = std::get<1>(inner_ball);
             }
-        #endif
-
+        }
         return _inner_ball;
     }
 
@@ -172,11 +203,18 @@ class HPolytope {
         return b;
     }
 
+    bool is_normalized ()
+    {
+        return normalized;
+    }
+
 
     // change the matrix A
     void set_mat(MT const& A2)
     {
         A = A2;
+        normalized = false;
+        has_ball = false;
     }
 
 
@@ -184,6 +222,7 @@ class HPolytope {
     void set_vec(VT const& b2)
     {
         b = b2;
+        has_ball = false;
     }
 
     Point get_mean_of_vertices() const
@@ -202,7 +241,7 @@ class HPolytope {
         std::cout << " " << A.rows() << " " << _d << " double" << std::endl;
         for (unsigned int i = 0; i < A.rows(); i++) {
             for (unsigned int j = 0; j < _d; j++) {
-                std::cout << A(i, j) << " ";
+                std::cout << A.coeff(i, j) << " ";
             }
             std::cout << "<= " << b(i) << std::endl;
         }
@@ -485,7 +524,7 @@ class HPolytope {
                                                      VT& Ar,
                                                      VT& Av,
                                                      NT const& lambda_prev,
-                                                     MT const& AA,
+                                                     DenseMT const& AA,
                                                      update_parameters& params) const
     {
 
@@ -501,7 +540,7 @@ class HPolytope {
         if(params.hit_ball) {
             Av.noalias() += (-2.0 * inner_prev) * (Ar / params.ball_inner_norm);
         } else {
-            Av.noalias() += (-2.0 * inner_prev) * AA.col(params.facet_prev);
+            Av.noalias() += ((-2.0 * inner_prev) * AA.col(params.facet_prev));
         }
         sum_nom.noalias() = b - Ar;
 
@@ -528,6 +567,56 @@ class HPolytope {
     }
 
 
+
+    template <typename update_parameters, typename set_type, typename AA_type>
+    std::pair<NT, int> line_positive_intersect(Point const& r,
+                                                     VT& Ar,
+                                                     VT& Av,
+                                                     NT const& lambda_prev,
+                                                     set_type& distances_set,
+                                                     AA_type const& AA,
+                                                     update_parameters& params) const
+    {
+        NT inner_prev = params.inner_vi_ak;
+        NT* Av_data = Av.data();
+
+        // Av += (-2.0 * inner_prev) * AA.col(params.facet_prev)
+
+        for (Eigen::SparseMatrix<double>::InnerIterator it(AA, params.facet_prev); it; ++it) {
+
+
+            // val(row) = (b(row) - Ar(row)) / Av(row) + params.moved_dist
+            // (val(row) - params.moved_dist) = (b(row) - Ar(row)) / Av(row)
+            // (val(row) - params.moved_dist) * Av(row) = b(row) - Ar(row)
+
+            *(Av_data + it.row()) += (-2.0 * inner_prev) * it.value();
+
+            // b(row) - Ar(row) = (old_val(row) - params.moved_dist) * old_Av(row) 
+            // new_val(row) = (b(row) - Ar(row)                                ) / new_Av(row) + params.moved_dist 
+            // new_val(row) = ((old_val(row) - params.moved_dist) * old_Av(row)) / new_Av(row) + params.moved_dist
+            
+            // new_val(row) = (old_val(row) - params.moved_dist) * old_Av(row)                                   / new_Av(row) + params.moved_dist;
+            // new_val(row) = (old_val(row) - params.moved_dist) * (new_Av(row) + 2.0 * inner_prev * it.value()) / new_Av(row) + params.moved_dist;
+            // new_val(row) = (old_val(row) - params.moved_dist) * (1 + (2.0 * inner_prev * it.value()) / new_Av(row) ) + params.moved_dist;
+
+            // new_val(row) = old_val(row) + (old_val(row) - params.moved_dist) * 2.0 * inner_prev * it.value() / new_Av(row)
+
+            // val(row) += (val(row) - params.moved_dist) * 2.0 * inner_prev * it.value() / *(Av_data + it.row());
+
+            NT val = distances_set.get_val(it.row());
+            val += (val - params.moved_dist) * 2.0 * inner_prev * it.value() / *(Av_data + it.row());
+            distances_set.change_val(it.row(), val, params.moved_dist);
+        }
+
+        std::pair<NT, int> ans = distances_set.get_min();
+        ans.first -= params.moved_dist;
+        params.inner_vi_ak = *(Av_data + ans.second);
+        params.facet_prev = ans.second;
+
+        return ans;
+    }
+
+
     template <typename update_parameters>
     std::pair<NT, int> line_positive_intersect(Point const& r,
                                                Point const& v,
@@ -541,7 +630,6 @@ class HPolytope {
 
         NT lamda = 0;
         VT sum_nom;
-        unsigned int j;
         int m = num_of_hyperplanes(), facet;
 
         Ar.noalias() += lambda_prev*Av;
@@ -623,12 +711,12 @@ class HPolytope {
 
         int m = num_of_hyperplanes();
 
-        lamdas.noalias() += A.col(rand_coord_prev)
-                         * (r_prev[rand_coord_prev] - r[rand_coord_prev]);
+        lamdas.noalias() += (DenseMT)(A.col(rand_coord_prev)
+                         * (r_prev[rand_coord_prev] - r[rand_coord_prev]));
         NT* data = lamdas.data();
 
         for (int i = 0; i < m; i++) {
-            NT a = A(i, rand_coord);
+            NT a = A.coeff(i, rand_coord);
 
             if (a == NT(0)) {
                 //std::cout<<"div0"<<std::endl;
@@ -819,9 +907,16 @@ class HPolytope {
 
 
     // Apply linear transformation, of square matrix T^{-1}, in H-polytope P:= Ax<=b
-    void linear_transformIt(MT const& T)
+    template<typename T_type>
+    void linear_transformIt(T_type const& T)
     {
-        A = A * T;
+        if constexpr (std::is_same<MT, DenseMT>::value) {
+            A = A * T;
+        } else {
+            A = (A * T).sparseView();
+        }
+        normalized = false;
+        has_ball = false;
     }
 
 
@@ -830,6 +925,7 @@ class HPolytope {
     void shift(const VT &c)
     {
         b -= A*c;
+        has_ball = false;
     }
 
 
@@ -859,13 +955,17 @@ class HPolytope {
 
     void normalize()
     {
+        if(normalized)
+            return;
         NT row_norm;
-        for (int i = 0; i < num_of_hyperplanes(); ++i)
-        {
+        for (int i = 0; i < A.rows(); ++i) {
             row_norm = A.row(i).norm();
-            A.row(i) = A.row(i) / row_norm;
-            b(i) = b(i) / row_norm;
+            if (row_norm != 0.0) {
+                A.row(i) /= row_norm;
+                b(i) /= row_norm;
+            }
         }
+        normalized = true;
     }
 
     void compute_reflection(Point& v, Point const&, int const& facet) const
@@ -903,10 +1003,34 @@ class HPolytope {
     }
 
     template <typename update_parameters>
-    void compute_reflection(Point &v, const Point &, update_parameters const& params) const {
+    void compute_reflection(Point &v, Point const&, update_parameters const& params) const {
+            Point a((-2.0 * params.inner_vi_ak) * A.row(params.facet_prev));
+            v += a;
+    }
+
+    // updates the velocity vector v and the position vector p after a reflection
+    // the real value of p is given by p + moved_dist * v
+    template <typename update_parameters>
+    auto compute_reflection(Point &v, Point &p, update_parameters const& params) const
+         -> std::enable_if_t<std::is_same_v<MT, Eigen::SparseMatrix<NT, Eigen::RowMajor>> && !std::is_same_v<update_parameters, int>, void> { // MT must be in RowMajor format
+            NT* v_data = v.pointerToData();
+            NT* p_data = p.pointerToData();
+            for(Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(A, params.facet_prev); it; ++it) {
+                *(v_data + it.col()) += (-2.0 * params.inner_vi_ak) * it.value();
+                *(p_data + it.col()) -= (-2.0 * params.inner_vi_ak * params.moved_dist) * it.value();
+            }
+    }
+
+    template <typename update_parameters>
+    NT compute_reflection(Point &v, const Point &, DenseMT 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&){}
diff --git a/src/volesti/include/convex_bodies/orderpolytope.h b/src/volesti/include/convex_bodies/orderpolytope.h
index 79f41e16..73b4e80f 100644
--- a/src/volesti/include/convex_bodies/orderpolytope.h
+++ b/src/volesti/include/convex_bodies/orderpolytope.h
@@ -94,12 +94,22 @@ class OrderPolytope {
         return _A.sparseView();
     }
 
+    // return the matrix A
+    MT get_dense_mat() const
+    {
+        return _A;
+    }
+
 
     VT get_vec() const
     {
         return b;
     }
 
+    bool is_normalized ()
+    {
+        return _normalized;
+    }
 
     // print polytope in Ax <= b format
     void print() const
diff --git a/src/volesti/include/convex_bodies/spectrahedra/spectrahedron.h b/src/volesti/include/convex_bodies/spectrahedra/spectrahedron.h
index e8f34138..37b997ea 100644
--- a/src/volesti/include/convex_bodies/spectrahedra/spectrahedron.h
+++ b/src/volesti/include/convex_bodies/spectrahedra/spectrahedron.h
@@ -341,6 +341,13 @@ class Spectrahedron {
     {
         return MT::Identity(lmi.dimension(), lmi.dimension());
     }
+    
+    bool is_normalized ()
+    {
+        return true;
+    }
+
+    void normalize() {}
 
     void resetFlags()
     {
diff --git a/src/volesti/include/convex_bodies/vpolyintersectvpoly.h b/src/volesti/include/convex_bodies/vpolyintersectvpoly.h
index a6dd97ca..ec2a7c5d 100644
--- a/src/volesti/include/convex_bodies/vpolyintersectvpoly.h
+++ b/src/volesti/include/convex_bodies/vpolyintersectvpoly.h
@@ -92,6 +92,10 @@ class IntersectionOfVpoly {
         return P1.get_vec();
     }
 
+    bool is_normalized () {
+        return true;
+    }
+
     MT get_T() const {
         return P1.get_mat();
     }
diff --git a/src/volesti/include/convex_bodies/vpolytope.h b/src/volesti/include/convex_bodies/vpolytope.h
index 82ca6689..9100756e 100644
--- a/src/volesti/include/convex_bodies/vpolytope.h
+++ b/src/volesti/include/convex_bodies/vpolytope.h
@@ -213,6 +213,10 @@ class VPolytope {
         return b;
     }
 
+    bool is_normalized () {
+        return true;
+    }
+
     // change the matrix V
     void set_mat(const MT &V2) {
         V = V2;
diff --git a/src/volesti/include/convex_bodies/zonoIntersecthpoly.h b/src/volesti/include/convex_bodies/zonoIntersecthpoly.h
index bbad43a6..e356db8d 100644
--- a/src/volesti/include/convex_bodies/zonoIntersecthpoly.h
+++ b/src/volesti/include/convex_bodies/zonoIntersecthpoly.h
@@ -65,6 +65,12 @@ class ZonoIntersectHPoly {
         return HP.get_vec();
     }
 
+    bool is_normalized () {
+        return true;
+    }
+
+    void normalize() {}
+
     std::pair<PointType,NT> InnerBall() const
     {
         return Z.InnerBall();
diff --git a/src/volesti/include/convex_bodies/zpolytope.h b/src/volesti/include/convex_bodies/zpolytope.h
index 1ad8868e..68ce0cf0 100644
--- a/src/volesti/include/convex_bodies/zpolytope.h
+++ b/src/volesti/include/convex_bodies/zpolytope.h
@@ -299,6 +299,11 @@ class Zonotope {
         return b;
     }
 
+    bool is_normalized ()
+    {
+        return true;
+    }
+
     // change the matrix V
     void set_mat(MT const& V2)
     {
diff --git a/src/volesti/include/generators/h_polytopes_generator.h b/src/volesti/include/generators/h_polytopes_generator.h
index 0c9293ba..aec3f1e8 100644
--- a/src/volesti/include/generators/h_polytopes_generator.h
+++ b/src/volesti/include/generators/h_polytopes_generator.h
@@ -10,6 +10,7 @@
 
 #include <exception>
 #include <chrono>
+#include <Eigen/Eigen>
 #include <boost/random/normal_distribution.hpp>
 #include <boost/random/uniform_real_distribution.hpp>
 
@@ -25,9 +26,9 @@
 template <class Polytope, class RNGType>
 Polytope random_hpoly(unsigned int dim, unsigned int m, int seed = std::numeric_limits<int>::signaling_NaN()) {
 
-    typedef typename Polytope::MT    MT;
     typedef typename Polytope::VT    VT;
     typedef typename Polytope::NT    NT;
+    typedef typename Eigen::Matrix<NT, Eigen::Dynamic, Eigen::Dynamic>    MT;
     typedef typename Polytope::PointType Point;
 
     int rng_seed = std::chrono::system_clock::now().time_since_epoch().count();
@@ -104,7 +105,7 @@ template <class Polytope, typename NT, class RNGType>
 Polytope skinny_random_hpoly(unsigned int dim, unsigned int m, const bool pre_rounding = false,
                              const NT eig_ratio = NT(1000.0), int seed = std::numeric_limits<int>::signaling_NaN()) {
 
-    typedef typename Polytope::MT    MT;
+    typedef typename Eigen::Matrix<NT, Eigen::Dynamic, Eigen::Dynamic>    MT;
     typedef typename Polytope::VT    VT;
     typedef typename Polytope::PointType Point;
 
diff --git a/src/volesti/include/generators/known_polytope_generators.h b/src/volesti/include/generators/known_polytope_generators.h
index 8d014621..e88963de 100644
--- a/src/volesti/include/generators/known_polytope_generators.h
+++ b/src/volesti/include/generators/known_polytope_generators.h
@@ -19,7 +19,7 @@
 template <class Polytope>
 Polytope generate_cube(const unsigned int& dim, const bool& Vpoly,typename Polytope::NT  scale=1) {
 
-    typedef typename Polytope::MT    MT;
+    typedef typename Eigen::Matrix<typename Polytope::NT, Eigen::Dynamic, Eigen::Dynamic>    MT;
     typedef typename Polytope::VT    VT;
     MT A;
     VT b;
@@ -84,7 +84,7 @@ template <typename Polytope>
 Polytope generate_cross(const unsigned int &dim, const bool &Vpoly) {
 
     unsigned int m;
-    typedef typename Polytope::MT    MT;
+    typedef typename Eigen::Matrix<typename Polytope::NT, Eigen::Dynamic, Eigen::Dynamic>    MT;
     typedef typename Polytope::VT    VT;
 
     MT A;
@@ -145,7 +145,7 @@ Polytope generate_cross(const unsigned int &dim, const bool &Vpoly) {
 /// @tparam Polytope Type of returned polytope
 template <typename Polytope>
 Polytope generate_simplex(const unsigned int &dim, const bool &Vpoly){
-    typedef typename Polytope::MT    MT;
+    typedef typename Eigen::Matrix<typename Polytope::NT, Eigen::Dynamic, Eigen::Dynamic>    MT;
     typedef typename Polytope::VT    VT;
 
     MT A;
@@ -198,7 +198,7 @@ Polytope generate_prod_simplex(const unsigned int &dim, bool Vpoly = false){
         return Perr;
     }
 
-    typedef typename Polytope::MT    MT;
+    typedef typename Eigen::Matrix<typename Polytope::NT, Eigen::Dynamic, Eigen::Dynamic>    MT;
     typedef typename Polytope::VT    VT;
 
     MT A;
@@ -274,7 +274,7 @@ Polytope generate_skinny_cube(const unsigned int &dim, bool Vpoly = false) {
         return Perr;
     }
 
-    typedef typename Polytope::MT    MT;
+    typedef typename Eigen::Matrix<typename Polytope::NT, Eigen::Dynamic, Eigen::Dynamic>    MT;
     typedef typename Polytope::VT    VT;
 
     MT A;
@@ -322,7 +322,7 @@ Polytope generate_birkhoff(unsigned int const& n) {
     unsigned int m = n * n;
     unsigned int d = n * n - 2 * n + 1;
 
-    typedef typename Polytope::MT MT;
+    typedef typename Eigen::Matrix<typename Polytope::NT, Eigen::Dynamic, Eigen::Dynamic> MT;
     typedef typename Polytope::VT VT;
 
     MT A = MT::Zero(m, d);
diff --git a/src/volesti/include/generators/order_polytope_generator.h b/src/volesti/include/generators/order_polytope_generator.h
index bfece556..739169e7 100644
--- a/src/volesti/include/generators/order_polytope_generator.h
+++ b/src/volesti/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_dense_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/src/volesti/include/lp_oracles/solve_lp.h b/src/volesti/include/lp_oracles/solve_lp.h
index f7013e99..a132ca99 100644
--- a/src/volesti/include/lp_oracles/solve_lp.h
+++ b/src/volesti/include/lp_oracles/solve_lp.h
@@ -80,8 +80,8 @@ std::pair<Point,NT> ComputeChebychevBall(const MT &A, const VT &b){
         sum=NT(0);
         for(j=0; j<d; j++){
             colno[j] = j+1;
-            row[j] = A(i,j);
-            sum+=A(i,j)*A(i,j);
+            row[j] = A.coeff(i,j);
+            sum+=A.coeff(i,j)*A.coeff(i,j);
         }
         colno[d] = d+1; /* last column */
         row[d] = std::sqrt(sum);
diff --git a/src/volesti/include/misc/poset.h b/src/volesti/include/misc/poset.h
index 69ba02ec..084bf66a 100644
--- a/src/volesti/include/misc/poset.h
+++ b/src/volesti/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/src/volesti/include/ode_solvers/implicit_midpoint.hpp b/src/volesti/include/ode_solvers/implicit_midpoint.hpp
index 786734d8..7f13bfe9 100644
--- a/src/volesti/include/ode_solvers/implicit_midpoint.hpp
+++ b/src/volesti/include/ode_solvers/implicit_midpoint.hpp
@@ -172,6 +172,8 @@ struct ImplicitMidpointODESolver {
       stream << '\n';
     }
   }
+
+  NT get_eta() const { return eta; }
 };
 
 #endif
diff --git a/src/volesti/include/ode_solvers/oracle_functors.hpp b/src/volesti/include/ode_solvers/oracle_functors.hpp
index 35fc8887..ab4546b6 100644
--- a/src/volesti/include/ode_solvers/oracle_functors.hpp
+++ b/src/volesti/include/ode_solvers/oracle_functors.hpp
@@ -394,4 +394,77 @@ struct HessianFunctor {
 
 };
 
+struct NonSphericalGaussianFunctor {
+  template <typename NT, typename Point>
+  struct parameters {
+    Point x0;
+    NT a;
+    NT eta;
+    unsigned int order;
+    NT L;
+    NT m;
+    NT kappa;
+    Eigen::Matrix<NT, Eigen::Dynamic, Eigen::Dynamic> inv_covariance_matrix;
+    
+    parameters(Point x0_, NT a_, NT eta_, Eigen::Matrix<NT, Eigen::Dynamic, Eigen::Dynamic> inv_covariance_matrix_)
+        : x0(x0_), a(a_), eta(eta_), order(2),
+          L(a_), m(a_), kappa(1),
+          inv_covariance_matrix(inv_covariance_matrix_) {};
+  };
+
+  template <typename Point>
+  struct GradientFunctor {
+    typedef typename Point::FT NT;
+    typedef std::vector<Point> pts;
+
+    parameters<NT, Point>& params;
+
+    GradientFunctor(parameters<NT, Point>& params_) : params(params_) {};
+
+    Point operator()(unsigned int const& i, pts const& xs, NT const& t) const {
+      if (i == params.order - 1) {
+          Point y = xs[0] - params.x0;
+          Eigen::Matrix<NT, Eigen::Dynamic, 1> result = -2.0 * params.a * params.inv_covariance_matrix * y.getCoefficients();
+          return Point(result);
+      } else {
+          return xs[i + 1];
+      }
+    }
+    Point operator()(Point const& x) {
+      Point y = x - params.x0;
+      Eigen::Matrix<NT, Eigen::Dynamic, 1> result = -2.0 * params.a * params.inv_covariance_matrix * y.getCoefficients();
+      return Point(result);
+    }
+  };
+
+  template <typename Point>
+  struct FunctionFunctor {
+    typedef typename Point::FT NT;
+
+    parameters<NT, Point>& params;
+
+    FunctionFunctor(parameters<NT, Point>& params_) : params(params_) {};
+
+    NT operator()(Point const& x) const {
+      Point y = x - params.x0;
+      Eigen::Matrix<NT, Eigen::Dynamic, 1> y_coeffs = y.getCoefficients();
+      return params.a * y_coeffs.dot(params.inv_covariance_matrix * y_coeffs);
+    }
+  };
+
+  template <typename Point>
+  struct HessianFunctor {
+    typedef typename Point::FT NT;
+
+    parameters<NT, Point>& params;
+
+    HessianFunctor(parameters<NT, Point>& params_) : params(params_) {};
+
+    Point operator()(Point const& x) const {
+        Eigen::Matrix<NT, Eigen::Dynamic, Eigen::Dynamic> result =   2.0 * params.a *  params.inv_covariance_matrix;
+        return Point(result);
+    }
+  };
+};
+
 #endif
diff --git a/src/volesti/include/preprocess/analytic_center_linear_ineq.h b/src/volesti/include/preprocess/analytic_center_linear_ineq.h
deleted file mode 100644
index 80d731f6..00000000
--- a/src/volesti/include/preprocess/analytic_center_linear_ineq.h
+++ /dev/null
@@ -1,134 +0,0 @@
-// VolEsti (volume computation and sampling library)
-
-// Copyright (c) 2024 Vissarion Fisikopoulos
-// Copyright (c) 2024 Apostolos Chalkis
-// Copyright (c) 2024 Elias Tsigaridas
-
-// Licensed under GNU LGPL.3, see LICENCE file
-
-
-#ifndef ANALYTIC_CENTER_H
-#define ANALYTIC_CENTER_H
-
-#include <tuple>
-
-#include "preprocess/max_inscribed_ball.hpp"
-#include "preprocess/feasible_point.hpp"
-#include "preprocess/mat_computational_operators.h"
-
-template <typename VT, typename NT>
-NT get_max_step(VT const& Ad, VT const& b_Ax)
-{
-    const int m = Ad.size();
-    NT max_element = std::numeric_limits<NT>::lowest(), max_element_temp;
-    for (int i = 0; i < m; i++) {
-        max_element_temp = Ad.coeff(i) / b_Ax.coeff(i);
-        if (max_element_temp > max_element) {
-            max_element = max_element_temp;
-        }
-    }
-
-    return NT(1) / max_element;
-}
-
-template <typename MT, typename VT, typename NT>
-void get_hessian_grad_logbarrier(MT const& A, MT const& A_trans, VT const& b, 
-                                 VT const& x, VT const& Ax, MT &H, VT &grad, VT &b_Ax)
-{
-    b_Ax.noalias() = b - Ax;
-    VT s = b_Ax.cwiseInverse();
-    VT s_sq = s.cwiseProduct(s);
-    // Gradient of the log-barrier function
-    grad.noalias() = A_trans * s;
-    // Hessian of the log-barrier function
-    update_Atrans_Diag_A<NT>(H, A_trans, A, s_sq.asDiagonal());
-}
-
-/*
-    This implementation computes the analytic center of a polytope given 
-    as a set of linear inequalities P = {x | Ax <= b}. The analytic center
-    is the minimizer of the log barrier function i.e., the optimal solution
-    of the following optimization problem (Convex Optimization, Boyd and Vandenberghe, Section 8.5.3),
-
-    \min -\sum \log(b_i - a_i^Tx), where a_i is the i-th row of A.
-    
-    The function solves the problem by using the Newton method.
-
-    Input: (i)   Matrix A, vector b such that the polytope P = {x | Ax<=b}
-           (ii)  The number of maximum iterations, max_iters
-           (iii) Tolerance parameter grad_err_tol to bound the L2-norm of the gradient
-           (iv)  Tolerance parameter rel_pos_err_tol to check the relative progress in each iteration
-
-    Output: (i)   The Hessian computed on the analytic center
-            (ii)  The analytic center of the polytope
-            (iii) A boolean variable that declares convergence
-    
-    Note: Using MT as to deal with both dense and sparse matrices, MT_dense will be the type of result matrix
-*/
-template <typename MT_dense, typename MT, typename VT, typename NT>
-std::tuple<MT_dense, VT, bool>  analytic_center_linear_ineq(MT const& A, VT const& b, VT const& x0,
-                                                            unsigned int const max_iters = 500,
-                                                            NT const grad_err_tol = 1e-08,
-                                                            NT const rel_pos_err_tol = 1e-12) 
-{
-    // Initialization
-    VT x = x0;
-    VT Ax = A * x;
-    const int n = A.cols(), m = A.rows();
-    MT H(n, n), A_trans = A.transpose();
-    VT grad(n), d(n), Ad(m), b_Ax(m), step_d(n), x_prev;
-    NT grad_err, rel_pos_err, rel_pos_err_temp, step;
-    unsigned int iter = 0;
-    bool converged = false;
-    const NT tol_bnd = NT(0.01);
-
-    auto llt = initialize_chol<NT>(A_trans, A);
-    get_hessian_grad_logbarrier<MT, VT, NT>(A, A_trans, b, x, Ax, H, grad, b_Ax);
-    
-    do {
-        iter++;
-        // Compute the direction
-        d.noalias() = - solve_vec<NT>(llt, H, grad);
-        Ad.noalias() = A * d;
-        // Compute the step length
-        step = std::min((NT(1) - tol_bnd) * get_max_step<VT, NT>(Ad, b_Ax), NT(1));
-        step_d.noalias() = step*d;
-        x_prev = x;
-        x += step_d;
-        Ax.noalias() += step*Ad;
-
-        // Compute the max_i\{ |step*d_i| ./ |x_i| \} 
-        rel_pos_err = std::numeric_limits<NT>::lowest();
-        for (int i = 0; i < n; i++)
-        {
-            rel_pos_err_temp = std::abs(step_d.coeff(i) / x_prev.coeff(i));
-            if (rel_pos_err_temp > rel_pos_err)
-            {
-                rel_pos_err = rel_pos_err_temp;
-            }
-        }
-        
-        get_hessian_grad_logbarrier<MT, VT, NT>(A, A_trans, b, x, Ax, H, grad, b_Ax);
-        grad_err = grad.norm();
-
-        if (iter >= max_iters || grad_err <= grad_err_tol || rel_pos_err <= rel_pos_err_tol)
-        {
-            converged = true;
-            break;
-        }
-    } while (true);
-    
-    return std::make_tuple(MT_dense(H), x, converged);
-}
-
-template <typename MT_dense, typename MT, typename VT, typename NT>
-std::tuple<MT_dense, VT, bool>  analytic_center_linear_ineq(MT const& A, VT const& b,
-                                                            unsigned int const max_iters = 500,
-                                                            NT const grad_err_tol = 1e-08,
-                                                            NT const rel_pos_err_tol = 1e-12) 
-{
-    VT x0 = compute_feasible_point(A, b);
-    return analytic_center_linear_ineq<MT_dense, MT, VT, NT>(A, b, x0, max_iters, grad_err_tol, rel_pos_err_tol);
-}
-
-#endif
diff --git a/src/volesti/include/preprocess/barrier_center_ellipsoid.hpp b/src/volesti/include/preprocess/barrier_center_ellipsoid.hpp
new file mode 100644
index 00000000..9f440f22
--- /dev/null
+++ b/src/volesti/include/preprocess/barrier_center_ellipsoid.hpp
@@ -0,0 +1,120 @@
+// VolEsti (volume computation and sampling library)
+
+// Copyright (c) 2024 Vissarion Fisikopoulos
+// Copyright (c) 2024 Apostolos Chalkis
+// Copyright (c) 2024 Elias Tsigaridas
+
+// Licensed under GNU LGPL.3, see LICENCE file
+
+
+#ifndef BARRIER_CENTER_ELLIPSOID_HPP
+#define BARRIER_CENTER_ELLIPSOID_HPP
+
+#include <tuple>
+
+#include "preprocess/max_inscribed_ball.hpp"
+#include "preprocess/feasible_point.hpp"
+#include "preprocess/rounding_util_functions.hpp"
+
+/*
+    This implementation computes the analytic or the volumetric center of a polytope given 
+    as a set of linear inequalities P = {x | Ax <= b}. The analytic center is the tminimizer
+    of the log barrier function, i.e., the optimal solution
+    of the following optimization problem (Convex Optimization, Boyd and Vandenberghe, Section 8.5.3),
+
+    \min -\sum \log(b_i - a_i^Tx), where a_i is the i-th row of A.
+    
+    The volumetric center is the minimizer of the volumetric barrier function, i.e., the optimal
+    solution of the following optimization problem,
+
+    \min logdet \nabla^2 f(x), where f(x) the log barrier function
+
+    The Vaidya center is the minimizer of the Vaidya barrier function, i.e., the optimal
+    solution of the following optimization problem,
+
+    \min logdet \nabla^2 f(x) + d/m f(x), where f(x) the log barrier function.
+    
+    The function solves the problems by using the Newton method.
+
+    Input: (i)   Matrix A, vector b such that the polytope P = {x | Ax<=b}
+           (ii)  The number of maximum iterations, max_iters
+           (iii) Tolerance parameter grad_err_tol to bound the L2-norm of the gradient
+           (iv)  Tolerance parameter rel_pos_err_tol to check the relative progress in each iteration
+
+    Output: (i)   The Hessian of the barrier function
+            (ii)  The analytic/volumetric center of the polytope
+            (iii) A boolean variable that declares convergence
+    
+    Note: Using MT as to deal with both dense and sparse matrices, MT_dense will be the type of result matrix
+*/
+template <typename MT_dense, int BarrierType, typename NT, typename MT, typename VT>
+std::tuple<MT_dense, VT, bool>  barrier_center_ellipsoid_linear_ineq(MT const& A, VT const& b, VT const& x0,
+                                                                     unsigned int const max_iters = 500,
+                                                                     NT const grad_err_tol = 1e-08,
+                                                                     NT const rel_pos_err_tol = 1e-12)
+{
+    // Initialization
+    VT x = x0;
+    VT Ax = A * x;
+    const int n = A.cols(), m = A.rows();
+    MT H(n, n), A_trans = A.transpose();
+    VT grad(n), d(n), Ad(m), b_Ax(m), step_d(n), x_prev;
+    NT grad_err, rel_pos_err, rel_pos_err_temp, step, obj_val, obj_val_prev;
+    unsigned int iter = 0;
+    bool converged = false;
+    const NT tol_bnd = NT(0.01), tol_obj = NT(1e-06);
+
+    auto [step_iter, max_step_multiplier] = init_step<BarrierType, NT>();
+    auto llt = initialize_chol<NT>(A_trans, A);
+    get_barrier_hessian_grad<MT_dense, BarrierType>(A, A_trans, b, x, Ax, llt,
+                                                    H, grad, b_Ax, obj_val);
+    do {
+        iter++;
+        // Compute the direction
+        d.noalias() = - solve_vec<NT>(llt, H, grad);
+        Ad.noalias() = A * d;
+        // Compute the step length
+        step = std::min(max_step_multiplier * get_max_step(Ad, b_Ax), step_iter);
+        step_d.noalias() = step*d;
+        x_prev = x;
+        x += step_d;
+        Ax.noalias() += step*Ad;
+
+        // Compute the max_i\{ |step*d_i| ./ |x_i| \} 
+        rel_pos_err = std::numeric_limits<NT>::lowest();
+        for (int i = 0; i < n; i++)
+        {
+            rel_pos_err_temp = std::abs(step_d.coeff(i) / x_prev.coeff(i));
+            if (rel_pos_err_temp > rel_pos_err)
+            {
+                rel_pos_err = rel_pos_err_temp;
+            }
+        }
+        
+        obj_val_prev = obj_val;
+        get_barrier_hessian_grad<MT_dense, BarrierType>(A, A_trans, b, x, Ax, llt,
+                                                        H, grad, b_Ax, obj_val);
+        grad_err = grad.norm();
+
+        if (iter >= max_iters || grad_err <= grad_err_tol || rel_pos_err <= rel_pos_err_tol)
+        {
+            converged = true;
+            break;
+        }
+        get_step_next_iteration<BarrierType>(obj_val_prev, obj_val, tol_obj, step_iter);
+    } while (true);
+    
+    return std::make_tuple(MT_dense(H), x, converged);
+}
+
+template <typename MT_dense, int BarrierType, typename NT, typename MT, typename VT>
+std::tuple<MT_dense, VT, bool>  barrier_center_ellipsoid_linear_ineq(MT const& A, VT const& b,
+                                                                     unsigned int const max_iters = 500,
+                                                                     NT const grad_err_tol = 1e-08,
+                                                                     NT const rel_pos_err_tol = 1e-12) 
+{
+    VT x0 = compute_feasible_point(A, b);
+    return barrier_center_ellipsoid_linear_ineq<MT_dense, BarrierType>(A, b, x0, max_iters, grad_err_tol, rel_pos_err_tol);
+}
+
+#endif // BARRIER_CENTER_ELLIPSOID_HPP
diff --git a/src/volesti/include/preprocess/crhmc/opts.h b/src/volesti/include/preprocess/crhmc/opts.h
index 573a69db..fffd29c5 100644
--- a/src/volesti/include/preprocess/crhmc/opts.h
+++ b/src/volesti/include/preprocess/crhmc/opts.h
@@ -44,6 +44,7 @@ template <typename Type> class opts {
   bool DynamicWeight = true; //Enable the use of dynamic weights for each variable when sampling
   bool DynamicStepSize = true;  // Enable adaptive step size that avoids low acceptance probability
   bool DynamicRegularizer = true; //Enable the addition of a regularization term
+  bool etaInitialize = true; // enable eta initialization
   Type regularization_factor=1e-20;
   /*Dynamic step choices*/
   Type warmUpStep = 10;
diff --git a/src/volesti/include/preprocess/estimate_L_smooth_parameter.hpp b/src/volesti/include/preprocess/estimate_L_smooth_parameter.hpp
index 180f2094..340ece75 100644
--- a/src/volesti/include/preprocess/estimate_L_smooth_parameter.hpp
+++ b/src/volesti/include/preprocess/estimate_L_smooth_parameter.hpp
@@ -21,7 +21,7 @@ template
     typename NegativeGradientFunctor,
     typename RandomNumberGenerator
 >
-double estimate_L_smooth(Polytope &P, Point &p, unsigned int const& walk_length, 
+double estimate_L_smooth(Polytope &P, Point &p, unsigned int const& walk_length,
                          NegativeGradientFunctor F, RandomNumberGenerator &rng)
 {
     typedef typename Point::FT NT;
@@ -44,11 +44,11 @@ double estimate_L_smooth(Polytope &P, Point &p, unsigned int const& walk_length,
     RandomWalk walk(P, p, rng);
     for (unsigned int i=0; i<rnum; ++i)
     {
-        walk.template apply(P, p, walk_length, rng);
+        walk.apply(P, p, walk_length, rng);
         randPoints[0] = p;
         listOfPoints.push_back(randPoints);
     }
-    
+
     NT L = std::numeric_limits<NT>::lowest(), Ltemp;
 
     for (auto pit=listOfPoints.begin(); pit!=(listOfPoints.end()-1); ++pit)
diff --git a/src/volesti/include/preprocess/inscribed_ellipsoid_rounding.hpp b/src/volesti/include/preprocess/inscribed_ellipsoid_rounding.hpp
index dd20d533..4b027f61 100644
--- a/src/volesti/include/preprocess/inscribed_ellipsoid_rounding.hpp
+++ b/src/volesti/include/preprocess/inscribed_ellipsoid_rounding.hpp
@@ -12,14 +12,9 @@
 #define INSCRIBED_ELLIPSOID_ROUNDING_HPP
 
 #include "preprocess/max_inscribed_ellipsoid.hpp"
-#include "preprocess/analytic_center_linear_ineq.h"
+#include "preprocess/barrier_center_ellipsoid.hpp"
 #include "preprocess/feasible_point.hpp"
 
-enum EllipsoidType
-{
-  MAX_ELLIPSOID = 1,
-  LOG_BARRIER = 2
-};
 
 template<typename MT, int ellipsoid_type, typename Custom_MT, typename VT, typename NT>
 inline static std::tuple<MT, VT, NT>
@@ -30,9 +25,11 @@ compute_inscribed_ellipsoid(Custom_MT A, VT b, VT const& x0,
     if constexpr (ellipsoid_type == EllipsoidType::MAX_ELLIPSOID)
     {
         return max_inscribed_ellipsoid<MT>(A, b, x0, maxiter, tol, reg);
-    } else if constexpr (ellipsoid_type == EllipsoidType::LOG_BARRIER)
+    } else if constexpr (ellipsoid_type == EllipsoidType::LOG_BARRIER ||
+                         ellipsoid_type == EllipsoidType::VOLUMETRIC_BARRIER ||
+                         ellipsoid_type == EllipsoidType::VAIDYA_BARRIER)
     {
-        return analytic_center_linear_ineq<MT, Custom_MT, VT, NT>(A, b, x0);
+        return barrier_center_ellipsoid_linear_ineq<MT, ellipsoid_type, NT>(A, b, x0);
     } else
     {
         std::runtime_error("Unknown rounding method.");
diff --git a/src/volesti/include/preprocess/max_inscribed_ball.hpp b/src/volesti/include/preprocess/max_inscribed_ball.hpp
index 573813dd..5f95970b 100644
--- a/src/volesti/include/preprocess/max_inscribed_ball.hpp
+++ b/src/volesti/include/preprocess/max_inscribed_ball.hpp
@@ -11,7 +11,7 @@
 #ifndef MAX_INSCRIBED_BALL_HPP
 #define MAX_INSCRIBED_BALL_HPP
 
-#include "preprocess/mat_computational_operators.h"
+#include "preprocess/rounding_util_functions.hpp"
 
 /*
     This implmentation computes the largest inscribed ball in a given convex polytope P.
diff --git a/src/volesti/include/preprocess/max_inscribed_ellipsoid.hpp b/src/volesti/include/preprocess/max_inscribed_ellipsoid.hpp
index f44f60ff..cc65006b 100644
--- a/src/volesti/include/preprocess/max_inscribed_ellipsoid.hpp
+++ b/src/volesti/include/preprocess/max_inscribed_ellipsoid.hpp
@@ -15,7 +15,7 @@
 
 #include <utility>
 #include <Eigen/Eigen>
-#include "preprocess/mat_computational_operators.h"
+#include "preprocess/rounding_util_functions.hpp"
 
 
 /*
diff --git a/src/volesti/include/preprocess/mat_computational_operators.h b/src/volesti/include/preprocess/rounding_util_functions.hpp
similarity index 65%
rename from src/volesti/include/preprocess/mat_computational_operators.h
rename to src/volesti/include/preprocess/rounding_util_functions.hpp
index af200d22..bb07d5b3 100644
--- a/src/volesti/include/preprocess/mat_computational_operators.h
+++ b/src/volesti/include/preprocess/rounding_util_functions.hpp
@@ -1,14 +1,15 @@
 // VolEsti (volume computation and sampling library)
 
-// Copyright (c) 2024 Vissarion Fisikopoulos
-// Copyright (c) 2024 Apostolos Chalkis
-// Copyright (c) 2024 Elias Tsigaridas
+// Copyright (c) 2012-2024 Vissarion Fisikopoulos
+// Copyright (c) 2018-2024 Apostolos Chalkis
+
+//Contributed and/or modified by Alexandros Manochis, as part of Google Summer of Code 2020 program.
 
 // Licensed under GNU LGPL.3, see LICENCE file
 
 
-#ifndef MAT_COMPUTATIONAL_OPERATORS_H
-#define MAT_COMPUTATIONAL_OPERATORS_H
+#ifndef ROUNDING_UTIL_FUNCTIONS_HPP
+#define ROUNDING_UTIL_FUNCTIONS_HPP
 
 #include <memory>
 
@@ -17,9 +18,36 @@
 #include "Spectra/include/Spectra/MatOp/SparseSymMatProd.h"
 
 
+enum EllipsoidType
+{
+  MAX_ELLIPSOID = 1,
+  LOG_BARRIER = 2,
+  VOLUMETRIC_BARRIER = 3,
+  VAIDYA_BARRIER = 4
+};
+
+template <int T>
+struct AssertBarrierFalseType : std::false_type {};
+
 template <typename T>
 struct AssertFalseType : std::false_type {};
 
+template <typename VT>
+auto get_max_step(VT const& Ad, VT const& b_Ax)
+{
+    using NT = typename VT::Scalar;
+    const int m = Ad.size();
+    NT max_element = std::numeric_limits<NT>::lowest(), max_element_temp;
+    for (int i = 0; i < m; i++) {
+        max_element_temp = Ad.coeff(i) / b_Ax.coeff(i);
+        if (max_element_temp > max_element) {
+            max_element = max_element_temp;
+        }
+    }
+
+    return NT(1) / max_element;
+}
+
 template <typename NT, typename MT>
 inline static auto
 initialize_chol(MT const& mat) 
@@ -29,7 +57,7 @@ initialize_chol(MT const& mat)
     if constexpr (std::is_same<MT, DenseMT>::value)
     {
         return std::make_unique<Eigen::LLT<MT>>();
-    } else if constexpr (std::is_same<MT, SparseMT>::value)  
+    } else if constexpr (std::is_base_of<Eigen::SparseMatrixBase<MT>, MT >::value)  
     {
         auto llt = std::make_unique<Eigen::SimplicialLLT<MT>>();
         llt->analyzePattern(mat);
@@ -50,7 +78,7 @@ initialize_chol(MT const& A_trans, MT const& A)
     if constexpr (std::is_same<MT, DenseMT>::value)
     {
         return std::make_unique<Eigen::LLT<MT>>();
-    } else if constexpr (std::is_same<MT, SparseMT>::value)  
+    } else if constexpr (std::is_base_of<Eigen::SparseMatrixBase<MT>, MT >::value)  
     {
         MT mat = A_trans * A;
         return initialize_chol<NT>(mat);
@@ -71,7 +99,7 @@ inline static VT solve_vec(std::unique_ptr<Eigen_lltMT> const& llt,
     {
         llt->compute(H);
         return llt->solve(b);
-    } else if constexpr (std::is_same<MT, SparseMT>::value)  
+    } else if constexpr (std::is_base_of<Eigen::SparseMatrixBase<MT>, MT >::value)  
     {
         llt->factorize(H);
         return llt->solve(b);
@@ -94,7 +122,7 @@ solve_mat(std::unique_ptr<Eigen_lltMT> const& llt,
         llt->compute(H);
         logdetE = llt->matrixL().toDenseMatrix().diagonal().array().log().sum();
         return llt->solve(mat);
-    } else if constexpr (std::is_same<MT, SparseMT>::value)  
+    } else if constexpr (std::is_base_of<Eigen::SparseMatrixBase<MT>, MT >::value)  
     {
         llt->factorize(H);
         logdetE = llt->matrixL().nestedExpression().diagonal().array().log().sum();
@@ -115,7 +143,7 @@ inline static void update_Atrans_Diag_A(MT &H, MT const& A_trans,
     if constexpr (std::is_same<MT, DenseMT>::value)
     {
         H.noalias() = A_trans * D * A;
-    } else if constexpr (std::is_same<MT, SparseMT>::value)  
+    } else if constexpr (std::is_base_of<Eigen::SparseMatrixBase<MT>, MT >::value)  
     {
         H = A_trans * D * A;
     } else 
@@ -133,7 +161,7 @@ inline static void update_Diag_A(MT &H, diag_MT const& D, MT const& A)
     if constexpr (std::is_same<MT, DenseMT>::value)
     {
         H.noalias() = D * A;
-    } else if constexpr (std::is_same<MT, SparseMT>::value)  
+    } else if constexpr (std::is_base_of<Eigen::SparseMatrixBase<MT>, MT >::value)  
     {
         H = D * A;
     } else 
@@ -151,7 +179,7 @@ inline static void update_A_Diag(MT &H, MT const& A, diag_MT const& D)
     if constexpr (std::is_same<MT, DenseMT>::value)
     {
         H.noalias() = A * D;
-    } else if constexpr (std::is_same<MT, SparseMT>::value)  
+    } else if constexpr (std::is_base_of<Eigen::SparseMatrixBase<MT>, MT >::value)  
     {
         H = A * D;
     } else 
@@ -170,7 +198,7 @@ get_mat_prod_op(MT const& E)
     if constexpr (std::is_same<MT, DenseMT>::value)
     {
         return std::make_unique<Spectra::DenseSymMatProd<NT>>(E);
-    } else if constexpr (std::is_same<MT, SparseMT>::value)  
+    } else if constexpr (std::is_base_of<Eigen::SparseMatrixBase<MT>, MT >::value)  
     {
         return std::make_unique<Spectra::SparseSymMatProd<NT>>(E);
     } else 
@@ -221,7 +249,7 @@ init_Bmat(MT &B, int const n, MT const& A_trans, MT const& A)
     if constexpr (std::is_same<MT, DenseMT>::value)
     {
         B.resize(n+1, n+1);
-    } else if constexpr (std::is_same<MT, SparseMT>::value)  
+    } else if constexpr (std::is_base_of<Eigen::SparseMatrixBase<MT>, MT >::value)  
     {
         // Initialize the structure of matrix B
         typedef Eigen::Triplet<NT> triplet;
@@ -271,7 +299,7 @@ update_Bmat(MT &B, VT const& AtDe, VT const& d,
         B.block(n, 0, 1, n).noalias() = AtDe.transpose();
         B(n, n) = d.sum();
         B.noalias() += 1e-14 * MT::Identity(n + 1, n + 1);
-    } else if constexpr (std::is_same<MT, SparseMT>::value)  
+    } else if constexpr (std::is_base_of<Eigen::SparseMatrixBase<MT>, MT >::value)  
     {
         MT AtD_A = AtD * A;
         int k = 0;
@@ -312,5 +340,91 @@ update_Bmat(MT &B, VT const& AtDe, VT const& d,
     }
 }
 
+template <int BarrierType, typename NT>
+std::tuple<NT, NT> init_step()
+{
+  if constexpr (BarrierType == EllipsoidType::LOG_BARRIER)
+  {
+    return {NT(1), NT(0.99)};
+  } else if constexpr (BarrierType == EllipsoidType::VOLUMETRIC_BARRIER ||
+                       BarrierType == EllipsoidType::VAIDYA_BARRIER)
+  {
+    return {NT(0.5), NT(0.4)};
+  } else {
+    static_assert(AssertBarrierFalseType<BarrierType>::value,
+            "Barrier type is not supported.");
+  }
+}
+
+template <typename MT_dense, int BarrierType, typename MT, typename VT, typename llt_type, typename NT>
+void get_barrier_hessian_grad(MT const& A, MT const& A_trans, VT const& b, 
+                              VT const& x, VT const& Ax, llt_type const& llt,
+                              MT &H, VT &grad, VT &b_Ax, NT &obj_val)
+{
+  b_Ax.noalias() = b - Ax;
+  VT s = b_Ax.cwiseInverse();
+  VT s_sq = s.cwiseProduct(s);
+  VT sigma;
+  // Hessian of the log-barrier function
+  update_Atrans_Diag_A<NT>(H, A_trans, A, s_sq.asDiagonal());
+
+  if constexpr (BarrierType == EllipsoidType::VOLUMETRIC_BARRIER ||
+                BarrierType == EllipsoidType::VAIDYA_BARRIER)
+  {
+    // Computing sigma(x)_i = (a_i^T H^{-1} a_i) / (b_i - a_i^Tx)^2
+    MT_dense HA = solve_mat(llt, H, A_trans, obj_val);
+    MT_dense aiHai = HA.transpose().cwiseProduct(A);
+    sigma = (aiHai.rowwise().sum()).cwiseProduct(s_sq);
+  }
+
+  if constexpr (BarrierType == EllipsoidType::LOG_BARRIER)
+  {
+    grad.noalias() = A_trans * s;
+  } else if constexpr (BarrierType == EllipsoidType::VOLUMETRIC_BARRIER)
+  {
+    // Gradient of the volumetric barrier function
+    grad.noalias() = A_trans * (s.cwiseProduct(sigma));
+    // Hessian of the volumetric barrier function
+    update_Atrans_Diag_A<NT>(H, A_trans, A, s_sq.cwiseProduct(sigma).asDiagonal());
+  } else if constexpr (BarrierType == EllipsoidType::VAIDYA_BARRIER)
+  {
+    const int m = b.size(), d = x.size();
+    NT const d_m = NT(d) / NT(m);
+    // Weighted gradient of the log barrier function
+    grad.noalias() = A_trans * s;
+    grad *= d_m;
+    // Add the gradient of the volumetric function
+    grad.noalias() += A_trans * (s.cwiseProduct(sigma));
+    // Weighted Hessian of the log barrier function
+    H *= d_m;
+    // Add the Hessian of the volumetric function
+    MT Hvol(d, d);
+    update_Atrans_Diag_A<NT>(Hvol, A_trans, A, s_sq.cwiseProduct(sigma).asDiagonal());
+    H += Hvol;
+    obj_val -= s.array().log().sum();
+  } else {
+    static_assert(AssertBarrierFalseType<BarrierType>::value,
+            "Barrier type is not supported.");
+  }
+}
+
+template <int BarrierType, typename NT>
+void get_step_next_iteration(NT const obj_val_prev, NT const obj_val,
+                             NT const tol_obj, NT &step_iter)
+{
+  if constexpr (BarrierType == EllipsoidType::LOG_BARRIER)
+  {
+    step_iter = NT(1);
+  } else if constexpr (BarrierType == EllipsoidType::VOLUMETRIC_BARRIER)
+  {
+    step_iter *= (obj_val_prev <= obj_val - tol_obj) ? NT(0.9) : NT(0.999);
+  } else if constexpr (BarrierType == EllipsoidType::VAIDYA_BARRIER)
+  {
+    step_iter *= NT(0.999);
+  } else {
+    static_assert(AssertBarrierFalseType<BarrierType>::value,
+            "Barrier type is not supported.");
+  }
+}
 
-#endif // MAT_COMPUTATIONAL_OPERATORS_H
+#endif // ROUNDING_UTIL_FUNCTIONS_HPP
\ No newline at end of file
diff --git a/src/volesti/include/random_walks/compute_diameter.hpp b/src/volesti/include/random_walks/compute_diameter.hpp
index b20aa839..71007d9e 100644
--- a/src/volesti/include/random_walks/compute_diameter.hpp
+++ b/src/volesti/include/random_walks/compute_diameter.hpp
@@ -30,11 +30,11 @@ struct compute_diameter
 };
 
 
-template <typename Point>
-struct compute_diameter<HPolytope<Point>>
+template <typename Point, typename MT>
+struct compute_diameter<HPolytope<Point, MT>>
 {
     template <typename NT>
-    static NT compute(HPolytope<Point> &P)
+    static NT compute(HPolytope<Point, MT> &P)
     {
         return NT(2) * std::sqrt(NT(P.dimension())) * P.InnerBall().second;
     }
diff --git a/src/volesti/include/random_walks/crhmc/additional_units/dynamic_step_size.hpp b/src/volesti/include/random_walks/crhmc/additional_units/dynamic_step_size.hpp
index 36721196..3daee1c3 100644
--- a/src/volesti/include/random_walks/crhmc/additional_units/dynamic_step_size.hpp
+++ b/src/volesti/include/random_walks/crhmc/additional_units/dynamic_step_size.hpp
@@ -50,7 +50,7 @@ class dynamic_step_size {
     consecutiveBadStep = IVT::Zero(simdLen);
     rejectSinceShrink = VT::Zero(simdLen);
 
-    if (options.warmUpStep > 0) {
+    if (options.warmUpStep > 0 && options.etaInitialize) {
       eta = 1e-3;
     } else {
       warmupFinished = true;
diff --git a/src/volesti/include/random_walks/crhmc/crhmc_walk.hpp b/src/volesti/include/random_walks/crhmc/crhmc_walk.hpp
index 21ed8d6b..9f1dd24a 100644
--- a/src/volesti/include/random_walks/crhmc/crhmc_walk.hpp
+++ b/src/volesti/include/random_walks/crhmc/crhmc_walk.hpp
@@ -164,6 +164,9 @@ struct CRHMCWalk {
     inline void disable_adaptive(){
       update_modules=false;
     }
+    NT get_current_eta() const {
+        return solver->get_eta();
+    }
     inline void apply(RandomNumberGenerator &rng,
       int walk_length = 1,
       bool metropolis_filter = true)
diff --git a/src/volesti/include/random_walks/gaussian_accelerated_billiard_walk.hpp b/src/volesti/include/random_walks/gaussian_accelerated_billiard_walk.hpp
index fbf43f44..532d7672 100644
--- a/src/volesti/include/random_walks/gaussian_accelerated_billiard_walk.hpp
+++ b/src/volesti/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,8 +24,6 @@
 
 #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
 {
@@ -57,56 +59,86 @@ struct GaussianAcceleratedBilliardWalk
 
 
     template
-            <
-                    typename Polytope,
-                    typename RandomNumberGenerator
-            >
+    <
+        typename Polytope,
+        typename RandomNumberGenerator,
+        typename E_type = typename Polytope::DenseMT
+    >
     struct Walk
     {
         typedef typename Polytope::PointType Point;
-        // typedef typename Polytope::MT MT;
+        typedef typename Polytope::DenseMT DenseMT;
+        typedef typename Polytope::VT VT;
         typedef typename Point::FT NT;
 
-        template <typename GenericPolytope, typename Ellipsoid>
+        void computeLcov(const E_type E)
+        {
+            if constexpr (std::is_base_of<Eigen::SparseMatrixBase<E_type>, E_type >::value) {
+                Eigen::SimplicialLLT<E_type> lltofE;
+                lltofE.compute(E);
+                if (lltofE.info() != Eigen::Success) {
+                    throw std::runtime_error("First Cholesky decomposition failed for sparse matrix!");
+                }
+                Eigen::SparseMatrix<NT> I(E.cols(), E.cols());
+                I.setIdentity();
+                Eigen::SparseMatrix<NT> E_inv = lltofE.solve(I);
+                Eigen::SimplicialLLT<E_type> lltofEinv;
+                lltofEinv.compute(E_inv);
+                if (lltofE.info() != Eigen::Success) {
+                    throw std::runtime_error("Second Cholesky decomposition failed for sparse matrix!");
+                }
+                _L_cov = lltofEinv.matrixL();
+            } else {
+                Eigen::LLT<E_type> lltOfE(E.llt().solve(E_type::Identity(E.cols(), E.cols()))); // compute the Cholesky decomposition of inv(E)
+                if (lltOfE.info() != Eigen::Success) {
+                    throw std::runtime_error("Cholesky decomposition failed for dense matrix!");
+                }
+                _L_cov = lltOfE.matrixL();
+            }
+        }
+
+        template <typename GenericPolytope>
         Walk(GenericPolytope& P,
              Point const& p,
-             Ellipsoid const& E,   // ellipsoid representing the Gaussian distribution
+             E_type 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);
+            computeLcov(E);
+            _E = E;
+            _AA.noalias() = (DenseMT)(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
+             E_type 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);
+            computeLcov(E);
+            _E = E;
+            _AA.noalias() = (DenseMT)(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 +146,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 +198,48 @@ 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() = (DenseMT)(P.get_mat() * _E);
+            _AEA = _AE.cwiseProduct((DenseMT)P.get_mat()).rowwise().sum();
+            /*
+            _AEA.resize(P.num_of_hyperplanes());
+            for(int i = 0; i < P.num_of_hyperplanes(); ++i)
+            {
+                _AEA(i) = _AE.row(i).dot(P.get_mat().row(i));
+            }*/
 
-            NT T = -std::log(rng.sample_urdist()) * _Len;
+            _v = GetDirection<Point>::apply(n, rng, false);
+            _v = Point(_L_cov.template triangularView<Eigen::Lower>() * _v.getCoefficients());
+
+            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 +251,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 +273,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
+        DenseMT _AA;
+        E_type _L_cov;   // LL' = inv(E)
+        DenseMT _AE;
+        E_type _E;
+        VT _AEA;
+        unsigned int _rho;
         update_parameters _update_parameters;
         typename Point::Coeff _lambdas;
         typename Point::Coeff _Av;
+        bool was_reset;
     };
 
 };
diff --git a/src/volesti/include/random_walks/gaussian_helpers.hpp b/src/volesti/include/random_walks/gaussian_helpers.hpp
index d1ab5429..76a15113 100644
--- a/src/volesti/include/random_walks/gaussian_helpers.hpp
+++ b/src/volesti/include/random_walks/gaussian_helpers.hpp
@@ -9,7 +9,14 @@ NT eval_exp(Point const& p, NT const& a)
 {
     return std::exp(-a * p.squared_length());
 }
-
+template <typename Point, typename NT, typename MT>
+NT eval_exp(Point const& p, MT const& inv_covariance_matrix, NT const& a_next, NT const& a_curr)
+{
+    Eigen::Matrix<NT, Eigen::Dynamic, 1> dist_vector = p.getCoefficients();
+    NT mahalanobis_dist = dist_vector.transpose() * inv_covariance_matrix * dist_vector;
+    NT log_ratio = (a_curr - a_next) * mahalanobis_dist;
+    return std::exp(log_ratio);
+}
 
 template <typename Point, typename NT>
 NT get_max(Point const& l, Point const& u, NT const& a_i)
diff --git a/src/volesti/include/random_walks/uniform_accelerated_billiard_walk.hpp b/src/volesti/include/random_walks/uniform_accelerated_billiard_walk.hpp
index 980c8e5d..f91abee4 100644
--- a/src/volesti/include/random_walks/uniform_accelerated_billiard_walk.hpp
+++ b/src/volesti/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
 
@@ -11,6 +12,126 @@
 #define RANDOM_WALKS_ACCELERATED_IMPROVED_BILLIARD_WALK_HPP
 
 #include "sampling/sphere.hpp"
+#include <Eigen/Eigen>
+#include <set>
+#include <vector>
+
+const double eps = 1e-10;
+
+// data structure which maintains the values of (b - Ar)/Av, and can extract the minimum positive value and the facet associated with it
+// vec[i].first contains the value of (b(i) - Ar(i))/Av(i) + moved_dist, where moved_dist is the total distance that the point has travelled so far
+// The heap will only contain the values from vec which are greater than moved_dist (so that they are positive)
+template<typename NT>
+class BoundaryOracleHeap {
+public:
+    int n, heap_size;
+    std::vector<std::pair<NT, int>> heap;
+    std::vector<std::pair<NT, int>> vec;
+
+private:
+    int siftDown(int index) {
+        while((index << 1) + 1 < heap_size) {
+            int child = (index << 1) + 1;
+            if(child + 1 < heap_size && heap[child + 1].first < heap[child].first - eps) {
+                child += 1;
+            }
+            if(heap[child].first < heap[index].first - eps)
+            {
+                std::swap(heap[child], heap[index]);
+                std::swap(vec[heap[child].second].second, vec[heap[index].second].second);
+                index = child;
+            } else {
+                return index;
+            }
+        }
+        return index;
+    }
+
+    int siftUp(int index) {
+        while(index > 0 && heap[(index - 1) >> 1].first - eps > heap[index].first) {
+            std::swap(heap[(index - 1) >> 1], heap[index]);
+            std::swap(vec[heap[(index - 1) >> 1].second].second, vec[heap[index].second].second);
+            index = (index - 1) >> 1;
+        }
+        return index;
+    }
+
+    // takes the index of a facet, and (in case it is in the heap) removes it from the heap.
+    void remove (int index) {
+        index = vec[index].second;
+        if(index == -1) {
+            return;
+        }
+        std::swap(heap[heap_size - 1], heap[index]);
+        std::swap(vec[heap[heap_size - 1].second].second, vec[heap[index].second].second);
+        vec[heap[heap_size - 1].second].second = -1;
+        heap_size -= 1;
+        index = siftDown(index);
+        siftUp(index);
+    }
+
+    // inserts a new value into the heap, with its associated facet
+    void insert (const std::pair<NT, int> val) {
+        vec[val.second].second = heap_size;
+        vec[val.second].first = val.first;
+        heap[heap_size++] = val;
+        siftUp(heap_size - 1);
+    }
+
+public:
+    BoundaryOracleHeap() {}
+
+    BoundaryOracleHeap(int n) : n(n), heap_size(0) {
+        heap.resize(n);
+        vec.resize(n);
+    }
+
+    // rebuilds the heap with the existing values from vec
+    // O(n)
+    void rebuild (const NT &moved_dist) {
+        heap_size = 0;
+        for(int i = 0; i < n; ++i) {
+            vec[i].second = -1;
+            if(vec[i].first - eps > moved_dist) {
+                vec[i].second = heap_size;
+                heap[heap_size++] = {vec[i].first, i};
+            }
+        }
+        for(int i = heap_size - 1; i >= 0; --i) {
+            siftDown(i);
+        }
+    }
+
+    // returns (b(i) - Ar(i))/Av(i) + moved_dist
+    // O(1)
+    NT get_val (const int &index) {
+        return vec[index].first;
+    }
+
+    // returns the nearest facet
+    // O(1)
+    std::pair<NT, int> get_min () {
+        return heap[0];
+    }
+
+    // changes the stored value for a given facet, and updates the heap accordingly
+    // O(logn)
+    void change_val(const int& index, const NT& new_val, const NT& moved_dist) {
+        if(new_val < moved_dist - eps) {
+            vec[index].first = new_val;
+            remove(index);
+        } else {
+            if(vec[index].second == -1) {
+                insert({new_val, index});
+            } else {
+                heap[vec[index].second].first = new_val;
+                vec[index].first = new_val;
+                siftDown(vec[index].second);
+                siftUp(vec[index].second);
+            }
+        }
+    }
+};
 
 
 // Billiard walk which accelarates each step for uniform distribution
@@ -37,35 +158,46 @@ struct AcceleratedBilliardWalk
     struct update_parameters
     {
         update_parameters()
-                :   facet_prev(0), hit_ball(false), inner_vi_ak(0.0), ball_inner_norm(0.0)
+                :   facet_prev(0), hit_ball(false), inner_vi_ak(0.0), ball_inner_norm(0.0), moved_dist(0.0)
         {}
         int facet_prev;
         bool hit_ball;
         double inner_vi_ak;
         double ball_inner_norm;
+        double moved_dist;
     };
 
     parameters param;
 
 
     template
-            <
-                    typename Polytope,
-                    typename RandomNumberGenerator
-            >
+    <
+            typename Polytope,
+            typename RandomNumberGenerator
+    >
     struct Walk
     {
         typedef typename Polytope::PointType Point;
         typedef typename Polytope::MT MT;
+        typedef typename Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> DenseMT;
         typedef typename Point::FT NT;
+        using AA_type = std::conditional_t< std::is_same_v<MT, typename Eigen::SparseMatrix<NT, Eigen::RowMajor>>, typename Eigen::SparseMatrix<NT>, DenseMT >; 
+        // AA is sparse colMajor if MT is sparse rowMajor, and Dense otherwise
 
         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);
-            _AA.noalias() = P.get_mat() * P.get_mat().transpose();
+            if constexpr (std::is_same<AA_type, Eigen::SparseMatrix<NT>>::value) {
+                _AA = (P.get_mat() * P.get_mat().transpose());
+            } else {
+                _AA.noalias() = (DenseMT)(P.get_mat() * P.get_mat().transpose());
+            }
             _rho = 1000 * P.dimension(); // upper bound for the number of reflections (experimental)
             initialize(P, p, rng);
         }
@@ -74,11 +206,18 @@ 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>
                                 ::template compute<NT>(P);
-            _AA.noalias() = P.get_mat() * P.get_mat().transpose();
+            if constexpr (std::is_same<AA_type, Eigen::SparseMatrix<NT>>::value) {
+                _AA = (P.get_mat() * P.get_mat().transpose());
+            } else {
+                _AA.noalias() = (DenseMT)(P.get_mat() * P.get_mat().transpose());
+            }
             _rho = 1000 * P.dimension(); // upper bound for the number of reflections (experimental)
             initialize(P, p, rng);
         }
@@ -96,6 +235,12 @@ struct AcceleratedBilliardWalk
             NT T;
             const NT dl = 0.995;
             int it;
+            typename Point::Coeff b;
+            NT* b_data;
+            if constexpr (std::is_same<MT, Eigen::SparseMatrix<NT, Eigen::RowMajor>>::value) {
+                b = P.get_vec();
+                b_data = b.data();
+            }
 
             for (auto j=0u; j<walk_length; ++j)
             {
@@ -104,8 +249,8 @@ 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 = P.line_first_positive_intersect(_p, _v, _lambdas, _Av, _update_parameters);
+                
                 if (T <= pbpair.first) {
                     _p += (T * _v);
                     _lambda_prev = T;
@@ -113,28 +258,52 @@ struct AcceleratedBilliardWalk
                 }
 
                 _lambda_prev = dl * pbpair.first;
-                _p += (_lambda_prev * _v);
+                if constexpr (std::is_same<MT, Eigen::SparseMatrix<NT, Eigen::RowMajor>>::value) {
+                    _update_parameters.moved_dist = _lambda_prev;
+                    NT* Ar_data = _lambdas.data();
+                    NT* Av_data = _Av.data();
+                    for(int i = 0; i < P.num_of_hyperplanes(); ++i) {
+                        _distances_set.vec[i].first = ( *(b_data + i) - (*(Ar_data + i)) ) / (*(Av_data + i));
+                    }
+                    // rebuild the heap with the new values of (b - Ar) / Av
+                    _distances_set.rebuild(_update_parameters.moved_dist);
+                } else {
+                    _p += (_lambda_prev * _v);
+                }
                 T -= _lambda_prev;
                 P.compute_reflection(_v, _p, _update_parameters);
                 it++;
 
                 while (it < _rho)
-                {
-                    std::pair<NT, int> pbpair
-                            = P.line_positive_intersect(_p, _v, _lambdas, _Av, _lambda_prev,
-                                                        _AA, _update_parameters);
+                {   
+                    std::pair<NT, int> pbpair;
+                    if constexpr (std::is_same<MT, Eigen::SparseMatrix<NT, Eigen::RowMajor>>::value) {
+                        pbpair = P.line_positive_intersect(_p, _lambdas, _Av, _lambda_prev,
+                                                           _distances_set, _AA, _update_parameters);
+                    } else {
+                        pbpair = P.line_positive_intersect(_p, _v, _lambdas, _Av, _lambda_prev,
+                                                           _AA, _update_parameters);
+                    }
                     if (T <= pbpair.first) {
                         _p += (T * _v);
                         _lambda_prev = T;
                         break;
                     }
                     _lambda_prev = dl * pbpair.first;
-                    _p += (_lambda_prev * _v);
+                    if constexpr (std::is_same<MT, Eigen::SparseMatrix<NT, Eigen::RowMajor>>::value) {
+                        _update_parameters.moved_dist += _lambda_prev;
+                    } else {
+                        _p += (_lambda_prev * _v);
+                    }
                     T -= _lambda_prev;
                     P.compute_reflection(_v, _p, _update_parameters);
                     it++;
                 }
-                if (it == _rho) _p = p0;
+                _p += _update_parameters.moved_dist * _v;
+                _update_parameters.moved_dist = 0.0;
+                if (it == _rho) {
+                    _p = p0;
+                }
             }
             p = _p;
         }
@@ -185,8 +354,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 +364,7 @@ struct AcceleratedBilliardWalk
             }
 
             if (Lmax > _L) {
-                if (P.dimension() <= 2500) 
-                {
+                if (P.dimension() <= 2500) {
                     update_delta(Lmax);
                 }
                 else{
@@ -235,6 +402,7 @@ struct AcceleratedBilliardWalk
             _Av.setZero(P.num_of_hyperplanes());
             _p = p;
             _v = GetDirection<Point>::apply(n, rng);
+            _distances_set = BoundaryOracleHeap<NT>(P.num_of_hyperplanes());
 
             NT T = -std::log(rng.sample_urdist()) * _L;
             Point p0 = _p;
@@ -296,11 +464,12 @@ struct AcceleratedBilliardWalk
         Point _p;
         Point _v;
         NT _lambda_prev;
-        MT _AA;
+        AA_type _AA;
         unsigned int _rho;
         update_parameters _update_parameters;
         typename Point::Coeff _lambdas;
         typename Point::Coeff _Av;
+        BoundaryOracleHeap<NT> _distances_set;
     };
 
 };
diff --git a/src/volesti/include/random_walks/uniform_accelerated_billiard_walk_parallel.hpp b/src/volesti/include/random_walks/uniform_accelerated_billiard_walk_parallel.hpp
index 1c2630f9..7ec437b2 100644
--- a/src/volesti/include/random_walks/uniform_accelerated_billiard_walk_parallel.hpp
+++ b/src/volesti/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/src/volesti/include/sampling/mmcs.hpp b/src/volesti/include/sampling/mmcs.hpp
index 07529ad6..02791535 100644
--- a/src/volesti/include/sampling/mmcs.hpp
+++ b/src/volesti/include/sampling/mmcs.hpp
@@ -81,7 +81,7 @@ bool perform_mmcs_step(Polytope &P,
         walk.template get_starting_point(P, p, q, 10, rng);
         for (int i = 0; i < window; i++)
         {
-            walk.template apply(P, q, walk_length, rng);
+            walk.apply(P, q, walk_length, rng);
             winPoints.col(i) = q.getCoefficients();
         }
         estimator.update_estimator(winPoints);
diff --git a/src/volesti/include/sampling/parallel_mmcs.hpp b/src/volesti/include/sampling/parallel_mmcs.hpp
index a8308f4c..da22e8f6 100644
--- a/src/volesti/include/sampling/parallel_mmcs.hpp
+++ b/src/volesti/include/sampling/parallel_mmcs.hpp
@@ -126,7 +126,7 @@ bool perform_parallel_mmcs_step(Polytope &P,
             walk.template get_starting_point(P, p, thread_random_walk_parameters, 10, rng);
             for (int i = 0; i < window; i++)
             {
-                walk.template apply(P, thread_random_walk_parameters, walk_length, rng);
+                walk.apply(P, thread_random_walk_parameters, walk_length, rng);
                 winPoints_per_thread[thread_index].col(i) = thread_random_walk_parameters.p.getCoefficients();
             }
 
diff --git a/src/volesti/include/sampling/random_point_generators.hpp b/src/volesti/include/sampling/random_point_generators.hpp
index e70e1461..998d8daf 100644
--- a/src/volesti/include/sampling/random_point_generators.hpp
+++ b/src/volesti/include/sampling/random_point_generators.hpp
@@ -35,7 +35,7 @@ struct RandomPointGenerator
         Walk walk(P, p, rng, parameters);
         for (unsigned int i=0; i<rnum; ++i)
         {
-            walk.template apply(P, p, walk_length, rng);
+            walk.apply(P, p, walk_length, rng);
             policy.apply(randPoints, p);
         }
     }
@@ -59,7 +59,7 @@ struct RandomPointGenerator
         Walk walk(P, p, rng);
         for (unsigned int i=0; i<rnum; ++i)
         {
-            walk.template apply(P, p, walk_length, rng);
+            walk.apply(P, p, walk_length, rng);
             policy.apply(randPoints, p);
         }
     }
@@ -95,7 +95,7 @@ struct MultivariateGaussianRandomPointGenerator
         Walk walk(P, p, E, rng, parameters);
         for (unsigned int i=0; i<rnum; ++i)
         {
-            walk.template apply(P, p, E, walk_length, rng);
+            walk.apply(P, p, walk_length, rng);
             policy.apply(randPoints, p);
         }
     }
@@ -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.apply(P, p, walk_length, rng);
             policy.apply(randPoints, p);
         }
     }
@@ -155,7 +155,7 @@ struct GaussianRandomPointGenerator
         Walk walk(P, p, a_i, rng);
         for (unsigned int i=0; i<rnum; ++i)
         {
-            walk.template apply(P, p, a_i, walk_length, rng);
+            walk.apply(P, p, a_i, walk_length, rng);
             policy.apply(randPoints, p);
         }
     }
@@ -184,7 +184,7 @@ struct GaussianRandomPointGenerator
 
         for (unsigned int i=0; i<rnum; ++i)
         {
-            walk.template apply(P, p, a_i, walk_length, rng);
+            walk.apply(P, p, a_i, walk_length, rng);
             policy.apply(randPoints, p);
         }
     }
@@ -215,7 +215,7 @@ struct BoundaryRandomPointGenerator
         Point p1(P.dimension()), p2(P.dimension());
         for (unsigned int i=0; i<rnum; ++i)
         {
-            walk.template apply(P, p1, p2, walk_length, rng);
+            walk.apply(P, p1, p2, walk_length, rng);
             policy.apply(randPoints, p1);
             policy.apply(randPoints, p2);
         }
@@ -344,7 +344,7 @@ struct ExponentialRandomPointGenerator
         bool success;
         for (unsigned int i=0; i<rnum; ++i)
         {
-            success = walk.template apply(P, p, walk_length, rng);
+            success = walk.apply(P, p, walk_length, rng);
             if (!success) {
                 //return;
                 throw std::range_error("A generated point is outside polytope");
@@ -379,7 +379,7 @@ struct ExponentialRandomPointGenerator
 
         for (unsigned int i=0; i<rnum; ++i)
         {
-            success = walk.template apply(P, p, walk_length, rng);
+            success = walk.apply(P, p, walk_length, rng);
             if (!success) {
                 //return;
                 throw std::range_error("A generated point is outside polytope");
diff --git a/src/volesti/include/sampling/random_point_generators_multithread.hpp b/src/volesti/include/sampling/random_point_generators_multithread.hpp
index 3e04a9b3..a4872fd7 100644
--- a/src/volesti/include/sampling/random_point_generators_multithread.hpp
+++ b/src/volesti/include/sampling/random_point_generators_multithread.hpp
@@ -19,7 +19,7 @@ template <typename GenericWalk>
 struct policy_storing
 {
     template <typename WalkPolicy, typename PointList, typename ThreadParameters>
-    static void store(WalkPolicy &policy, PointList &randPoints, ThreadParameters &thread_random_walk_parameters) 
+    static void store(WalkPolicy &policy, PointList &randPoints, ThreadParameters &thread_random_walk_parameters)
     {
         policy.apply(randPoints, thread_random_walk_parameters.p);
     }
@@ -30,7 +30,7 @@ template <>
 struct policy_storing<BRDHRWalk_multithread>
 {
     template <typename WalkPolicy, typename PointList, typename ThreadParameters>
-    static void store(WalkPolicy &policy, PointList &randPoints, ThreadParameters &thread_random_walk_parameters) 
+    static void store(WalkPolicy &policy, PointList &randPoints, ThreadParameters &thread_random_walk_parameters)
     {
         policy.apply(randPoints, thread_random_walk_parameters.p1);
         policy.apply(randPoints, thread_random_walk_parameters.p2);
@@ -38,7 +38,7 @@ struct policy_storing<BRDHRWalk_multithread>
 };
 
 template <>
-struct policy_storing<BCDHRWalk_multithread> : policy_storing<BRDHRWalk_multithread> 
+struct policy_storing<BCDHRWalk_multithread> : policy_storing<BRDHRWalk_multithread>
 {};
 
 
@@ -72,7 +72,7 @@ struct RandomPointGeneratorMultiThread
 
         omp_set_num_threads(num_threads);
         unsigned int d = P.dimension(), m = P.num_of_hyperplanes();
-        
+
         std::vector<int> num_points_per_thread(rnum%num_threads, rnum/num_threads+1);
         std::vector<int> a(num_threads - rnum%num_threads, rnum/num_threads);
         num_points_per_thread.insert(num_points_per_thread.end(), a.begin(), a.end());
@@ -88,7 +88,7 @@ struct RandomPointGeneratorMultiThread
 
             for (unsigned int it = 0; it < num_points_per_thread[thread_index]; it++)
             {
-                walk.template apply(P, thread_random_walk_parameters, walk_length, rng);
+                walk.apply(P, thread_random_walk_parameters, walk_length, rng);
                 #pragma omp critical
                 {
                     policy_storing<Walk>::template store(policy, randPoints, thread_random_walk_parameters);
@@ -135,7 +135,7 @@ struct RandomPointGeneratorMultiThread
 
             for (unsigned int it = 0; it < num_points_per_thread[thread_index]; it++)
             {
-                walk.template apply(P, thread_random_walk_parameters, walk_length, rng);
+                walk.apply(P, thread_random_walk_parameters, walk_length, rng);
                 #pragma omp critical
                 {
                     policy_storing<Walk>::template store(policy, randPoints, thread_random_walk_parameters);
@@ -178,7 +178,7 @@ struct GaussianPointGeneratorMultiThread
 
         omp_set_num_threads(num_threads);
         unsigned int d = P.dimension(), m = P.num_of_hyperplanes();
-        
+
         std::vector<int> num_points_per_thread(rnum%num_threads, rnum/num_threads+1);
         std::vector<int> a(num_threads - rnum%num_threads, rnum/num_threads);
         num_points_per_thread.insert(num_points_per_thread.end(), a.begin(), a.end());
@@ -194,7 +194,7 @@ struct GaussianPointGeneratorMultiThread
 
             for (unsigned int it = 0; it < num_points_per_thread[thread_index]; it++)
             {
-                walk.template apply(P, thread_random_walk_parameters, a_i, walk_length, rng);
+                walk.apply(P, thread_random_walk_parameters, a_i, walk_length, rng);
                 #pragma omp critical
                 {
                     policy_storing<Walk>::template store(policy, randPoints, thread_random_walk_parameters);
@@ -226,7 +226,7 @@ struct GaussianPointGeneratorMultiThread
 
         omp_set_num_threads(num_threads);
         unsigned int d = P.dimension(), m = P.num_of_hyperplanes();
-        
+
         std::vector<int> num_points_per_thread(rnum%num_threads, rnum/num_threads+1);
         std::vector<int> a(num_threads - rnum%num_threads, rnum/num_threads);
         num_points_per_thread.insert(num_points_per_thread.end(), a.begin(), a.end());
@@ -242,7 +242,7 @@ struct GaussianPointGeneratorMultiThread
 
             for (unsigned int it = 0; it < num_points_per_thread[thread_index]; it++)
             {
-                walk.template apply(P, thread_random_walk_parameters, a_i, walk_length, rng);
+                walk.apply(P, thread_random_walk_parameters, a_i, walk_length, rng);
                 #pragma omp critical
                 {
                     policy_storing<Walk>::template store(policy, randPoints, thread_random_walk_parameters);
diff --git a/src/volesti/include/sampling/sample_correlation_matrices.hpp b/src/volesti/include/sampling/sample_correlation_matrices.hpp
index 50b7a45d..b2bd118b 100644
--- a/src/volesti/include/sampling/sample_correlation_matrices.hpp
+++ b/src/volesti/include/sampling/sample_correlation_matrices.hpp
@@ -14,6 +14,21 @@
 
 #include <sampling/sampling.hpp>
 
+template<typename NT, typename MT>
+MT getCoefficientsFromMatrix(const MT& mat) {
+    int n = mat.rows();
+    int d = n * (n - 1) / 2;
+    Eigen::Matrix<NT, Eigen::Dynamic, 1> coeffs(d);
+    int k = 0;
+    for (int i = 0; i < n; ++i) {
+        for (int j = 0; j < i; ++j) {
+            coeffs(k) = mat(i, j);
+            ++k;
+    	}
+    }
+    return coeffs;
+}
+
 // New implementations for sampling correlation matrices with CorrelationSpectrahedron and CorrelationSpectrahedron_MT
 
 template
diff --git a/src/volesti/include/volume/volume_cooling_balls.hpp b/src/volesti/include/volume/volume_cooling_balls.hpp
index d90604fa..7f3d00d5 100644
--- a/src/volesti/include/volume/volume_cooling_balls.hpp
+++ b/src/volesti/include/volume/volume_cooling_balls.hpp
@@ -466,7 +466,7 @@ NT estimate_ratio(PolyBall1 &Pb1,
 
     do
     {
-        walk.template apply(Pb1, p, walk_length, rng);
+        walk.apply(Pb1, p, walk_length, rng);
     } while(!estimate_ratio_generic(Pb2, p, error, ratio_parameters));
 
     return NT(ratio_parameters.count_in) / NT(ratio_parameters.tot_count);
@@ -674,14 +674,14 @@ NT estimate_ratio_interval(PolyBall1 &Pb1,
 
     for (int i = 0; i < ratio_parameters.W; ++i)
     {
-        walk.template apply(Pb1, p, walk_length, rng);
+        walk.apply(Pb1, p, walk_length, rng);
         full_sliding_window(Pb2, p, ratio_parameters);
     }
 
     ratio_parameters.mean = ratio_parameters.sum / NT(ratio_parameters.W);
 
     do {
-        walk.template apply(Pb1, p, walk_length, rng);
+        walk.apply(Pb1, p, walk_length, rng);
     } while (!estimate_ratio_interval_generic(Pb2, p, error, zp, ratio_parameters));
 
     return NT(ratio_parameters.count_in) / NT(ratio_parameters.tot_count);
diff --git a/src/volesti/include/volume/volume_cooling_gaussians.hpp b/src/volesti/include/volume/volume_cooling_gaussians.hpp
index ca06ce31..62835bd8 100644
--- a/src/volesti/include/volume/volume_cooling_gaussians.hpp
+++ b/src/volesti/include/volume/volume_cooling_gaussians.hpp
@@ -241,7 +241,7 @@ void compute_annealing_schedule(Polytope& P,
         // Compute some ratios to decide if this is the last gaussian
         for (unsigned  int j = 0; j < steps; j++)
         {
-            walk.template apply(P, p, a_vals[it], walk_length, rng);
+            walk.apply(P, p, a_vals[it], walk_length, rng);
             curr_its += 1.0;
             curr_fn += eval_exp(p, next_a) / eval_exp(p, a_vals[it]);
         }
@@ -398,7 +398,7 @@ double volume_cooling_gaussians(Polytope& Pin,
 
         while (!done || (*itsIt)<min_steps)
         {
-            walk.template apply(P, p, *avalsIt, walk_length, rng);
+            walk.apply(P, p, *avalsIt, walk_length, rng);
 
             *itsIt = *itsIt + 1.0;
             *fnIt = *fnIt + eval_exp(p,*(avalsIt+1)) / eval_exp(p,*avalsIt);
diff --git a/src/volesti/include/volume/volume_cooling_gaussians_crhmc.hpp b/src/volesti/include/volume/volume_cooling_gaussians_crhmc.hpp
new file mode 100644
index 00000000..779f5b58
--- /dev/null
+++ b/src/volesti/include/volume/volume_cooling_gaussians_crhmc.hpp
@@ -0,0 +1,459 @@
+// VolEsti (volume computation and sampling library)
+
+// Copyright (c) 2012-2024 Vissarion Fisikopoulos
+// Copyright (c) 2018-2024 Apostolos Chalkis
+// Copyright (c) 2024 Vladimir Necula
+
+// Contributed and/or modified by Vladimir Necula, as part of Google Summer of
+// Code 2024 program.
+
+// Licensed under GNU LGPL.3, see LICENCE file
+
+// References
+// Yunbum Kook, Yin Tat Lee, Ruoqi Shen, Santosh S. Vempala. "Sampling with
+// Riemannian Hamiltonian
+// Monte Carlo in a Constrained Space"
+#ifndef VOLUME_COOLING_GAUSSIANS_CRHMC_HPP
+#define VOLUME_COOLING_GAUSSIANS_CRHMC_HPP
+
+//#define VOLESTI_DEBUG
+
+#include "volume/volume_cooling_gaussians.hpp"
+#include "preprocess/crhmc/crhmc_problem.h"
+
+template
+<
+    typename CRHMCWalkType,
+    typename crhmc_walk_params,
+    int simdLen,
+    typename Grad,
+    typename Func,
+    typename CrhmcProblem,
+    typename Polytope,
+    typename Point,
+    typename NT,
+    typename RandomNumberGenerator
+>
+NT get_next_gaussian(Polytope& P,
+                     Point &p,
+                     NT const& a,
+                     const unsigned int &N,
+                     const NT &ratio,
+                     const NT &C,
+                     const unsigned int& walk_length,
+                     RandomNumberGenerator& rng,
+                     Grad& g,
+                     Func& f,
+                     crhmc_walk_params& parameters,
+                     CrhmcProblem& problem,
+                     CRHMCWalkType& crhmc_walk)
+{
+
+    NT last_a = a;
+    NT last_ratio = 0.1;
+    //k is needed for the computation of the next variance a_{i+1} = a_i * (1-1/d)^k
+    NT k = 1.0;
+    const NT tol = 0.00001;
+    bool done=false;
+    std::vector<NT> fn(N,NT(0.0));
+    std::list<Point> randPoints;
+    typedef typename std::vector<NT>::iterator viterator;
+
+    //sample N points
+    PushBackWalkPolicy push_back_policy;
+    bool raw_output = false; 
+
+    typedef CrhmcRandomPointGenerator<CRHMCWalkType> CRHMCRandomPointGenerator;
+
+    CRHMCRandomPointGenerator::apply(problem, p, N, walk_length, randPoints,
+                                    push_back_policy, rng, g, f, parameters, crhmc_walk, simdLen, raw_output);
+
+    while (!done)
+    {
+        NT new_a = last_a * std::pow(ratio,k);
+
+        auto fnit = fn.begin();
+        for (auto pit=randPoints.begin(); pit!=randPoints.end(); ++pit, fnit++)
+        {
+            *fnit = eval_exp(*pit, new_a)/eval_exp(*pit, last_a);
+        }
+        std::pair<NT, NT> mv = get_mean_variance(fn);
+
+        // Compute a_{i+1}
+        if (mv.second/(mv.first * mv.first)>=C || mv.first/last_ratio<1.0+tol)
+        {
+            if (k != 1.0)
+            {
+                k = k / 2;
+            }
+            done = true;
+        } 
+        else {
+            k = 2 * k;
+        }
+        last_ratio = mv.first;
+    }
+    return last_a * std::pow(ratio, k);
+}
+
+// Compute the sequence of spherical gaussians
+template
+<
+    int simdLen,
+    typename Polytope,
+    typename NT,
+    typename RandomNumberGenerator
+>
+void compute_annealing_schedule(Polytope& P,
+                                NT const& ratio,
+                                NT const& C,
+                                NT const& frac,
+                                unsigned int const& N,
+                                unsigned int const& walk_length,
+                                NT const& chebychev_radius,
+                                NT const& error,
+                                std::vector<NT>& a_vals,
+                                RandomNumberGenerator& rng)
+{
+    typedef typename Polytope::PointType Point;
+    typedef typename Polytope::MT 	MT;
+    
+    typedef typename GaussianFunctor::FunctionFunctor<Point>    Func;
+    typedef typename GaussianFunctor::GradientFunctor<Point>    Grad;
+    typedef typename GaussianFunctor::HessianFunctor<Point>     Hess;
+    typedef typename GaussianFunctor::parameters<NT, Point>     func_params;
+
+    typedef crhmc_input<MT, Point, Func, Grad, Hess> Input;
+    typedef crhmc_problem<Point, Input> CrhmcProblem;   
+
+    typedef ImplicitMidpointODESolver<Point, NT, CrhmcProblem, Grad, simdLen> Solver;
+
+    typedef typename CRHMCWalk::template Walk
+            <
+                    Point,
+                    CrhmcProblem,
+                    RandomNumberGenerator,
+                    Grad,
+                    Func,
+                    Solver
+            > CRHMCWalkType;
+
+    typedef typename CRHMCWalk::template parameters
+          <
+                  NT,
+                  Grad
+          > crhmc_walk_params;
+
+    typedef CrhmcRandomPointGenerator<CRHMCWalkType> CRHMCRandomPointGenerator;
+
+
+    // Compute the first gaussian
+    // This uses the function from the standard volume_cooling_gaussians.hpp
+    get_first_gaussian(P, frac, chebychev_radius, error, a_vals);
+    NT a_stop = 0.0;
+    const NT tol = 0.001;
+    unsigned int it = 0;
+    unsigned int n = P.dimension();
+    const unsigned int totalSteps = ((int)150/((1.0 - frac) * error))+1;
+
+    if (a_vals[0]<a_stop) a_vals[0] = a_stop;
+
+#ifdef VOLESTI_DEBUG
+    std::cout<<"first gaussian computed "<< a_vals[0] <<std::endl;
+#endif
+
+#ifdef VOLESTI_DEBUG
+    std::cout<<"Computing the sequence of gaussians..\n"<<std::endl;
+#endif
+
+    while (true)
+    {
+
+        NT curr_fn = 0;
+        NT curr_its = 0;
+        auto steps = totalSteps;
+
+        //TODO: potential problem creation and preprocessing optimization
+        
+        //create the crhmc problem for this variance
+        int dimension = P.dimension();
+        func_params f_params = func_params(Point(dimension), a_vals[it], 1);
+        
+        Func f(f_params);
+        Grad g(f_params);
+        Hess h(f_params);
+
+        Input input = convert2crhmc_input<Input, Polytope, Func, Grad, Hess>(P, f, g, h);
+
+        typedef crhmc_problem<Point, Input> CrhmcProblem;
+        CrhmcProblem problem = CrhmcProblem(input);
+
+        Point p = Point(problem.center);
+
+        if(problem.terminate){return;}
+
+        problem.options.simdLen = simdLen;
+        crhmc_walk_params params(input.df, p.dimension(), problem.options);
+
+        if (input.df.params.eta > 0) {
+            params.eta = input.df.params.eta;
+        }
+
+        int dim;
+        dim = p.dimension();
+
+        //create the walk object for this problem
+        CRHMCWalkType walk = CRHMCWalkType(problem, p, input.df, input.f, params);
+
+        // Compute the next gaussian
+        NT next_a = get_next_gaussian<CRHMCWalkType, crhmc_walk_params, simdLen, Grad, Func, CrhmcProblem>
+                      (P, p, a_vals[it], N, ratio, C, walk_length, rng, g, f, params, problem, walk);
+
+#ifdef VOLESTI_DEBUG
+    std::cout<<"Next Gaussian " << next_a <<std::endl;
+#endif
+
+        // Compute some ratios to decide if this is the last gaussian
+        for (unsigned  int j = 0; j < steps; j++)
+        {
+            walk.apply(rng, walk_length);
+            p = walk.getPoint();
+#ifdef VOLESTI_DEBUG
+    std::cout<<"Walk point " << std::endl;
+    p.print();
+#endif
+            curr_its += 1.0;
+            curr_fn += eval_exp(p, next_a) / eval_exp(p, a_vals[it]);
+        }
+
+#ifdef VOLESTI_DEBUG
+    std::cout<<"Condition function ratio " << curr_fn/curr_its <<std::endl;
+#endif
+
+        // Remove the last gaussian.
+        // Set the last a_i equal to zero
+        if (next_a>0 && curr_fn/curr_its>(1.0+tol))
+        {
+            a_vals.push_back(next_a);
+            it++;
+        } else if (next_a <= 0)
+        {
+            a_vals.push_back(a_stop);
+            it++;
+            break;
+        } else {
+            a_vals[it] = a_stop;
+            break;
+        }
+    }
+#ifdef VOLESTI_DEBUG
+    std::cout<<"first gaussian after WHILE "<< a_vals[0] <<std::endl;
+#endif
+}
+
+template
+<
+    typename Polytope,
+    typename RandomNumberGenerator,
+    int simdLen = 8
+>
+double volume_cooling_gaussians(Polytope& Pin,
+                                RandomNumberGenerator& rng,
+                                double const& error = 0.1,
+                                unsigned int const& walk_length = 1)
+{
+    typedef typename Polytope::PointType Point;
+    typedef typename Point::FT 	NT;
+    typedef typename Polytope::VT 	VT;
+    typedef typename Polytope::MT 	MT;
+    typedef typename GaussianFunctor::FunctionFunctor<Point>    Func;
+    typedef typename GaussianFunctor::GradientFunctor<Point>    Grad;
+    typedef typename GaussianFunctor::HessianFunctor<Point>     Hess;
+    typedef typename GaussianFunctor::parameters<NT, Point>     func_params;
+
+    typedef crhmc_input<MT, Point, Func, Grad, Hess> Input;
+    typedef crhmc_problem<Point, Input> CrhmcProblem;   
+
+    typedef ImplicitMidpointODESolver<Point, NT, CrhmcProblem, Grad, simdLen> Solver;
+
+    typedef typename CRHMCWalk::template Walk
+            <
+                    Point,
+                    CrhmcProblem,
+                    RandomNumberGenerator,
+                    Grad,
+                    Func,
+                    Solver
+            > CRHMCWalkType;
+
+    typedef typename CRHMCWalk::template parameters
+          <
+                  NT,
+                  Grad
+          > crhmc_walk_params;
+    
+    typedef CrhmcRandomPointGenerator<CRHMCWalkType> CRHMCRandomPointGenerator;
+
+    auto P(Pin); //copy and work with P because we are going to shift
+    unsigned int n = P.dimension();
+    unsigned int m = P.num_of_hyperplanes();
+    gaussian_annealing_parameters<NT> parameters(P.dimension());
+
+    // Consider Chebychev center as an internal point
+    auto InnerBall = P.ComputeInnerBall();
+    if (InnerBall.second < 0.0) return -1.0;
+
+    Point c = InnerBall.first;
+    NT radius = InnerBall.second;
+
+    // Move the chebychev center to the origin and apply the same shifting to the polytope
+    P.shift(c.getCoefficients());
+
+    // Computing the sequence of gaussians
+#ifdef VOLESTI_DEBUG
+    std::cout<<"\n\nComputing annealing...\n"<<std::endl;
+    double tstart2 = (double)clock()/(double)CLOCKS_PER_SEC;
+#endif
+
+    // Initialization for the schedule annealing
+    std::vector<NT> a_vals;
+    NT ratio = parameters.ratio;
+    NT C = parameters.C;
+    unsigned int N = parameters.N;
+
+    compute_annealing_schedule<simdLen>(P, ratio, C, parameters.frac, N, walk_length, radius, error, a_vals, rng);
+
+#ifdef VOLESTI_DEBUG
+    std::cout<<"All the variances of schedule_annealing computed in = "
+            << (double)clock()/(double)CLOCKS_PER_SEC-tstart2<<" sec"<<std::endl;
+    auto j=0;
+    for (auto avalIt = a_vals.begin(); avalIt!=a_vals.end(); avalIt++, j++){
+        std::cout<<"a_"<<j<<" = "<<*avalIt<<" ";
+    }
+    std::cout<<std::endl<<std::endl;
+#endif
+
+    // Initialization for the approximation of the ratios
+    unsigned int W = parameters.W;
+    unsigned int mm = a_vals.size()-1;
+    std::vector<NT> last_W2(W,0);
+    std::vector<NT> fn(mm,0);
+    std::vector<NT> its(mm,0);
+    VT lamdas;
+    lamdas.setZero(m);
+    NT vol = std::pow(M_PI/a_vals[0], (NT(n))/2.0);
+    unsigned int i=0;
+
+    typedef typename std::vector<NT>::iterator viterator;
+    viterator itsIt = its.begin();
+    viterator avalsIt = a_vals.begin();
+    viterator minmaxIt;
+
+
+#ifdef VOLESTI_DEBUG
+    std::cout<<"volume of the first gaussian = "<<vol<<"\n"<<std::endl;
+    std::cout<<"computing ratios..\n"<<std::endl;
+#endif
+
+    //iterate over the number of ratios
+    for (viterator fnIt = fn.begin();
+         fnIt != fn.end();
+         fnIt++, itsIt++, avalsIt++, i++)
+    {
+        //initialize convergence test
+        bool done = false;
+        NT curr_eps = error/std::sqrt((NT(mm)));
+        NT min_val = std::numeric_limits<NT>::min();
+        NT max_val = std::numeric_limits<NT>::max();
+        unsigned int min_index = W-1;
+        unsigned int max_index = W-1;
+        unsigned int index = 0;
+        unsigned int min_steps = 0;
+        std::vector<NT> last_W = last_W2;
+
+        // Set the radius for the ball walk
+        //creating the walk object
+        int dimension = P.dimension();
+        func_params f_params = func_params(Point(dimension), *avalsIt, 1);
+        
+        Func f(f_params);
+        Grad g(f_params);
+        Hess h(f_params);
+
+        //create the crhmc problem
+        Input input = convert2crhmc_input<Input, Polytope, Func, Grad, Hess>(P, f, g, h);
+
+        CrhmcProblem problem = CrhmcProblem(input);
+
+        Point p = Point(problem.center);
+
+        if(problem.terminate){return 0;}
+
+        problem.options.simdLen=simdLen;
+        crhmc_walk_params params(input.df, p.dimension(), problem.options);
+
+        if (input.df.params.eta > 0) {
+            params.eta = input.df.params.eta;
+        }
+
+        CRHMCWalkType walk = CRHMCWalkType(problem, p, input.df, input.f, params);
+
+        while (!done || (*itsIt)<min_steps)
+        {
+            walk.apply(rng, walk_length);
+            p = walk.getPoint();
+            *itsIt = *itsIt + 1.0;
+            *fnIt = *fnIt + eval_exp(p,*(avalsIt+1)) / eval_exp(p,*avalsIt);
+            NT val = (*fnIt) / (*itsIt);
+
+            last_W[index] = val;
+            if (val <= min_val)
+            {
+                min_val = val;
+                min_index = index;
+            } else if (min_index == index)
+            {
+                minmaxIt = std::min_element(last_W.begin(), last_W.end());
+                min_val = *minmaxIt;
+                min_index = std::distance(last_W.begin(), minmaxIt);
+            }
+
+            if (val >= max_val)
+            {
+                max_val = val;
+                max_index = index;
+            } else if (max_index == index)
+            {
+                minmaxIt = std::max_element(last_W.begin(), last_W.end());
+                max_val = *minmaxIt;
+                max_index = std::distance(last_W.begin(), minmaxIt);
+            }
+
+            if ( (max_val-min_val)/max_val <= curr_eps/2.0 )
+            {
+                done=true;
+            }
+
+            index = index%W + 1;
+            if (index == W) index = 0;
+        }
+#ifdef VOLESTI_DEBUG
+        std::cout << "ratio " << i << " = " << (*fnIt) / (*itsIt)
+                  << " N_" << i << " = " << *itsIt << std::endl;
+#endif
+        vol *= ((*fnIt) / (*itsIt));
+    }
+
+#ifdef VOLESTI_DEBUG
+        NT sum_of_steps = 0.0;
+        for(viterator it = its.begin(); it != its.end(); ++it) {
+            sum_of_steps += *it;
+        }
+        auto steps= int(sum_of_steps);
+        std::cout<<"\nTotal number of steps = "<<steps<<"\n"<<std::endl;
+#endif
+
+    return vol;
+}
+
+#endif // VOLUME_COOLING_GAUSSIANS_CRHMC_HPP
diff --git a/src/volesti/include/volume/volume_cooling_nonspherical_gaussians_crhmc.hpp b/src/volesti/include/volume/volume_cooling_nonspherical_gaussians_crhmc.hpp
new file mode 100644
index 00000000..fd5ce185
--- /dev/null
+++ b/src/volesti/include/volume/volume_cooling_nonspherical_gaussians_crhmc.hpp
@@ -0,0 +1,594 @@
+// VolEsti (volume computation and sampling library)
+
+// Copyright (c) 2012-2024 Vissarion Fisikopoulos
+// Copyright (c) 2018-2024 Apostolos Chalkis
+// Copyright (c) 2024 Vladimir Necula
+
+// Contributed and/or modified by Vladimir Necula, as part of Google Summer of
+// Code 2024 program.
+
+// Licensed under GNU LGPL.3, see LICENCE file
+
+#ifndef VOLUME_COOLING_NON_SPHERICAL_GAUSSIANS_CRHMC_HPP
+#define VOLUME_COOLING_NON_SPHERICAL_GAUSSIANS_CRHMC_HPP
+
+//#define VOLESTI_DEBUG
+
+#include "volume/volume_cooling_gaussians.hpp"
+#include "preprocess/crhmc/crhmc_problem.h"
+#include "preprocess/max_inscribed_ellipsoid.hpp"
+#include "preprocess/inscribed_ellipsoid_rounding.hpp"
+
+////////////////////////////// Algorithms
+
+// Gaussian Anealling
+
+// Compute a_{i+1} when a_i is given
+
+template
+<
+    typename CRHMCWalkType,
+    typename crhmc_walk_params,
+    int simdLen,
+    typename Grad,
+    typename Func,
+    typename CrhmcProblem,
+    typename Point,
+    typename NT,
+    typename RandomNumberGenerator
+>
+void burnIn(Point &p,
+            const unsigned int& walk_length,
+            RandomNumberGenerator& rng,
+            Grad& g,
+            Func& f,
+            crhmc_walk_params& parameters,
+            CrhmcProblem& problem,
+            CRHMCWalkType& crhmc_walk,
+            NT burnIn_sample) 
+{
+    std::list<Point> randPoints;
+
+    // burnIn
+    PushBackWalkPolicy push_back_policy;
+    bool raw_output = false;
+    typedef CrhmcRandomPointGenerator<CRHMCWalkType> CRHMCRandomPointGenerator;
+
+    CRHMCRandomPointGenerator::apply(problem, p, burnIn_sample, walk_length, randPoints,
+                                    push_back_policy, rng, g, f, parameters, crhmc_walk, simdLen, raw_output);
+
+}
+
+
+template
+<
+    int simdLen,
+    typename Polytope,
+    typename Point,
+    typename NT,
+    typename MT,
+    typename RandomNumberGenerator
+>
+NT get_next_gaussian(Polytope& P,
+                    Point& p,
+                    NT const& a,
+                    const unsigned int &N,
+                    const NT &ratio,
+                    const NT &C,
+                    const unsigned int& walk_length,
+                    RandomNumberGenerator& rng,
+                    MT const& inv_covariance_matrix) 
+{
+    typedef typename NonSphericalGaussianFunctor::FunctionFunctor<Point>    Func;
+    typedef typename NonSphericalGaussianFunctor::GradientFunctor<Point>    Grad;
+    typedef typename NonSphericalGaussianFunctor::HessianFunctor<Point>     Hess;
+    typedef typename NonSphericalGaussianFunctor::parameters<NT, Point>     func_params;
+
+    typedef crhmc_input<MT, Point, Func, Grad, ZeroFunctor<Point>> Input;
+
+    typedef crhmc_problem<Point, Input> CrhmcProblem;   
+
+    typedef ImplicitMidpointODESolver<Point, NT, CrhmcProblem, Grad, simdLen> Solver;
+
+    typedef typename CRHMCWalk::template Walk
+            <
+                    Point,
+                    CrhmcProblem,
+                    RandomNumberGenerator,
+                    Grad,
+                    Func,
+                    Solver
+            > CRHMCWalkType;
+
+    typedef typename CRHMCWalk::template parameters
+          <
+                  NT,
+                  Grad
+          > crhmc_walk_params;
+
+    typedef CrhmcRandomPointGenerator<CRHMCWalkType> CRHMCRandomPointGenerator;
+
+    // Create the CRHMC problem for this variance
+    int dimension = P.dimension();
+    func_params f_params = func_params(Point(dimension), a, -1, inv_covariance_matrix);
+
+    Func f(f_params);
+    Grad g(f_params);
+    Hess h(f_params);
+    ZeroFunctor<Point> zerof;
+
+    Input input = convert2crhmc_input<Input, Polytope, Func, Grad, ZeroFunctor<Point>>(P, f, g, zerof);
+
+    typedef crhmc_problem<Point, Input> CrhmcProblem;
+    
+    CrhmcProblem problem = CrhmcProblem(input);
+
+    if(problem.terminate) { return 0;}
+
+    problem.options.simdLen = simdLen;
+
+    crhmc_walk_params params(input.df, p.dimension(), problem.options);
+
+    // Create the walk object for this problem
+    CRHMCWalkType walk = CRHMCWalkType(problem, p, input.df, input.f, params);
+
+    burnIn<CRHMCWalkType, crhmc_walk_params, simdLen, Grad, Func, CrhmcProblem>(
+        p, walk_length, rng, g, f, params, problem, walk, 2.5 * N);
+
+    NT last_a = a;
+    NT last_ratio = 0.1;
+    NT k = 1.0;
+    const NT tol = 0.00001;
+    bool done = false;
+    std::vector<NT> fn(N, NT(0.0));
+    std::list<Point> randPoints;
+
+    // sample N points
+    PushBackWalkPolicy push_back_policy;
+    bool raw_output = false;
+    typedef CrhmcRandomPointGenerator<CRHMCWalkType> CRHMCRandomPointGenerator;
+
+    CRHMCRandomPointGenerator::apply(problem, p, N, walk_length, randPoints,
+                                    push_back_policy, rng, g, f, params, walk, simdLen, raw_output);
+
+    while (!done) {
+        NT new_a = last_a * std::pow(ratio, k);
+
+        auto fnit = fn.begin();
+        for (auto pit = randPoints.begin(); pit != randPoints.end(); ++pit, fnit++) {
+            *fnit = eval_exp(*pit, inv_covariance_matrix, new_a, last_a);
+        }
+
+        std::pair<NT, NT> mv = get_mean_variance(fn);
+
+        // Compute a_{i+1}
+        if (mv.second / (mv.first * mv.first) >= C || mv.first / last_ratio < 1.0 + tol) {
+            if (k != 1.0) {
+                k = k / 2;
+            }
+            done = true;
+        } else {
+            k = 2 * k;
+        }
+
+        last_ratio = mv.first;
+    }
+
+    // Return the new a value as a scalar
+    return last_a * std::pow(ratio, k);
+}
+
+// Compute the sequence of non spherical gaussians
+template<
+    int simdLen,
+    typename Point,
+    typename Polytope,
+    typename NT,
+    typename MT,
+    typename RandomNumberGenerator
+>
+void compute_annealing_schedule(Polytope Pin_copy,
+                                Polytope P_copy,
+                                NT const& ratio,
+                                NT const& C,
+                                NT const& frac,
+                                unsigned int const& N,
+                                unsigned int const& walk_length,
+                                NT const& error,
+                                std::vector<NT>& a_vals,
+                                MT const& inv_covariance_matrix,
+                                RandomNumberGenerator& rng) 
+{
+    
+    typedef typename NonSphericalGaussianFunctor::FunctionFunctor<Point>    Func;
+    typedef typename NonSphericalGaussianFunctor::GradientFunctor<Point>    Grad;
+    typedef typename NonSphericalGaussianFunctor::HessianFunctor<Point>     Hess;
+    typedef typename NonSphericalGaussianFunctor::parameters<NT, Point>     func_params;
+
+    typedef crhmc_input<MT, Point, Func, Grad, ZeroFunctor<Point>> Input;
+
+    typedef crhmc_problem<Point, Input> CrhmcProblem;   
+
+    typedef ImplicitMidpointODESolver<Point, NT, CrhmcProblem, Grad, simdLen> Solver;
+
+    typedef typename CRHMCWalk::template Walk
+            <
+                    Point,
+                    CrhmcProblem,
+                    RandomNumberGenerator,
+                    Grad,
+                    Func,
+                    Solver
+            > CRHMCWalkType;
+
+    typedef typename CRHMCWalk::template parameters
+          <
+                  NT,
+                  Grad
+          > crhmc_walk_params;
+
+    typedef CrhmcRandomPointGenerator<CRHMCWalkType> CRHMCRandomPointGenerator;
+
+    // compute the first variance
+    // for this we need P
+    auto ball = P_copy.ComputeInnerBall();
+    P_copy.shift(ball.first.getCoefficients()); // when using max_ellipsoid for rounding this center is the origin, but if we use other covariances this is different than the origin
+    get_first_gaussian(P_copy, frac, ball.second, error, a_vals); // the function included from volume_cooling_gaussians.hpp (spherical gaussians)
+#ifdef VOLESTI_DEBUG
+    std::cout << "first gaussian computed " << a_vals[0] << std::endl;
+#endif
+
+    //for the rest we need Pin
+    NT a_stop = 0.0;
+    const NT tol = 0.001;
+    unsigned int it = 0;
+    unsigned int n = Pin_copy.dimension();
+    const unsigned int totalSteps = ((int)150/((1.0 - frac) * error)) + 1;
+
+    if (a_vals[0]<a_stop) a_vals[0] = a_stop;
+
+#ifdef VOLESTI_DEBUG
+    std::cout << "Computing the sequence of gaussians..\n" << std::endl;
+#endif
+    NT initial_eta;
+    Point start_point;
+
+    int dimension = Pin_copy.dimension();
+    func_params initial_f_params = func_params(Point(dimension), a_vals[0], -1, inv_covariance_matrix);
+    Func initial_f(initial_f_params);
+    Grad initial_g(initial_f_params);
+    ZeroFunctor<Point> initial_zerof;
+
+    Input initial_input = convert2crhmc_input<Input, Polytope, Func, Grad, ZeroFunctor<Point>>(Pin_copy, initial_f, initial_g, initial_zerof);
+    CrhmcProblem initial_problem = CrhmcProblem(initial_input);
+
+    Point initial_p = Point(initial_problem.center);
+    initial_problem.options.simdLen = simdLen;
+    crhmc_walk_params initial_params(initial_input.df, initial_p.dimension(), initial_problem.options);
+    CRHMCWalkType initial_walk = CRHMCWalkType(initial_problem, initial_p, initial_input.df, initial_input.f, initial_params);
+    
+    burnIn<CRHMCWalkType, crhmc_walk_params, simdLen, Grad, Func, CrhmcProblem>(
+        initial_p, walk_length, rng, initial_g, initial_f, initial_params, initial_problem, initial_walk, 2.5 * N);
+
+    //fix eta
+    initial_eta = initial_walk.get_current_eta();
+
+    //fix point
+    start_point = initial_problem.center;
+
+    while (true) {
+
+        // Compute the next gaussian
+        NT next_a = get_next_gaussian<simdLen>(
+            Pin_copy, start_point, a_vals[it], N, ratio, C, walk_length, rng, inv_covariance_matrix);
+
+        // Decide if this is the last one
+        NT curr_fn = 0;
+        NT curr_its = 0;
+        auto steps = totalSteps;
+
+        //TODO: potential problem creation and preprocessing optimization
+
+        // Create the CRHMC problem for this variance
+        int dimension = Pin_copy.dimension();
+        func_params f_params = func_params(Point(dimension), a_vals[it], -1, inv_covariance_matrix);
+
+        Func f(f_params);
+        Grad g(f_params);
+        ZeroFunctor<Point> zerof;
+
+        Input input = convert2crhmc_input<Input, Polytope, Func, Grad, ZeroFunctor<Point>>(Pin_copy, f, g, zerof);
+
+        typedef crhmc_problem<Point, Input> CrhmcProblem;
+        
+        CrhmcProblem problem = CrhmcProblem(input);
+
+       if(problem.terminate) { return; }
+
+        problem.options.simdLen = simdLen;
+
+        crhmc_walk_params params(input.df, start_point.dimension(), problem.options);
+
+        //eta initialization with the previous walk eta
+        problem.options.etaInitialize = false;
+        params.eta = initial_eta;
+
+        // Create the walk object for this problem
+        CRHMCWalkType walk = CRHMCWalkType(problem, start_point, input.df, input.f, params);
+
+        // Compute some ratios to decide if this is the last gaussian
+        for (unsigned int j = 0; j < steps; j++) 
+        {
+            walk.apply(rng, walk_length);
+            start_point = walk.getPoint();
+            curr_its += 1.0;    
+            curr_fn += eval_exp(start_point, inv_covariance_matrix, next_a, a_vals[it]);
+        }
+
+        //restore the new eta and start point, by looking at the walk after its operations
+        initial_eta = walk.get_current_eta();  
+
+        // Remove the last gaussian.
+        // Set the last a_i equal to zero
+        if (next_a > 0 && curr_fn / curr_its > (1.0 + tol)) 
+        {
+            a_vals.push_back(next_a);
+            it++;
+        } else if (next_a <= 0) 
+        {
+            a_vals.push_back(a_stop);
+            it++;
+            break;
+        } else {
+            a_vals[it] = a_stop;
+            break;
+        }
+    }
+}
+
+
+template <typename NT>
+struct non_gaussian_annealing_parameters
+{
+    non_gaussian_annealing_parameters(unsigned int d)
+        : frac(0.1)
+        , ratio(NT(1) - NT(1) / NT(d))
+        , C(NT(2))
+        , N(500 * ((int) C) + ((int) (d * d)))
+        , W(6 * d * d + 800)
+    {}
+
+    NT frac;
+    NT ratio;
+    NT C;
+    unsigned int N;
+    unsigned int W;
+};
+
+
+template
+<
+    typename Polytope,
+    typename RandomNumberGenerator,
+    typename WalkTypePolicy = CRHMCWalk,
+    int simdLen = 8
+>
+double non_spherical_crhmc_volume_cooling_gaussians(Polytope& Pin,
+                                RandomNumberGenerator& rng,
+                                double const& error = 0.1,
+                                unsigned int const& walk_length = 1)
+{
+    typedef typename Polytope::PointType Point;
+    typedef typename Point::FT     NT;
+    typedef typename Polytope::VT  VT;
+    typedef typename Polytope::MT  MT;
+    typedef typename NonSphericalGaussianFunctor::FunctionFunctor<Point>    Func;
+    typedef typename NonSphericalGaussianFunctor::GradientFunctor<Point>    Grad;
+    typedef typename NonSphericalGaussianFunctor::parameters<NT, Point>     func_params;
+
+    typedef crhmc_input<MT, Point, Func, Grad, ZeroFunctor<Point>> Input;
+    typedef crhmc_problem<Point, Input> CrhmcProblem;   
+
+    typedef ImplicitMidpointODESolver<Point, NT, CrhmcProblem, Grad, simdLen> Solver;
+
+    typedef typename WalkTypePolicy::template Walk
+            <
+                    Point,
+                    CrhmcProblem,
+                    RandomNumberGenerator,
+                    Grad,
+                    Func,
+                    Solver
+            > CRHMCWalkType;
+
+    typedef typename WalkTypePolicy::template parameters
+          <
+                  NT,
+                  Grad
+          > crhmc_walk_params;
+
+    typedef HPolytope<Point> HPOLYTOPE;
+
+    HPOLYTOPE P(Pin.dimension(), Pin.get_mat(), Pin.get_vec());
+    HPOLYTOPE newPin(Pin.dimension(), Pin.get_mat(), Pin.get_vec());
+    unsigned int n = P.dimension();
+    unsigned int m = P.num_of_hyperplanes();
+    
+    //compute inscribed ellipsoid
+    NT tol = std::pow(10, -6.0), reg = std::pow(10, -4.0);
+    unsigned int maxiter = 100;
+    P.normalize();
+    VT x0 = compute_feasible_point(P.get_mat(), P.get_vec());
+    auto ellipsoid_result = compute_inscribed_ellipsoid<MT, EllipsoidType::MAX_ELLIPSOID>(P.get_mat(), P.get_vec(), x0, maxiter, tol, reg);
+
+    // extract the covariance matrix and the center of the ellipsoid
+    MT inv_covariance_matrix = std::get<0>(ellipsoid_result); //this is the covariance to use in the telescopic product
+    VT center = std::get<1>(ellipsoid_result);
+    MT covariance_matrix = inv_covariance_matrix.inverse(); 
+
+    newPin.shift(center); //we shift the initial polytope so that the origin is the center of the gaussians
+    P.shift(center);
+
+    // we apply the rounding transformation
+    Eigen::LLT<MT> lltOfA(covariance_matrix);
+    auto L = lltOfA.matrixL();
+    P.linear_transformIt(L);
+
+
+    // Initialize the gaussian_annealing_parameters struct
+    non_gaussian_annealing_parameters<NT> parameters(P.dimension());
+
+    // Computing the sequence of gaussians
+#ifdef VOLESTI_DEBUG
+    std::cout<<"\n\nComputing annealing...\n"<<std::endl;
+    double tstart2 = (double)clock()/(double)CLOCKS_PER_SEC;
+#endif
+
+    // Initialization for the schedule annealing
+    std::vector<NT> a_vals;
+    NT ratio = parameters.ratio;
+    NT C = parameters.C;
+    unsigned int N = parameters.N;
+
+    compute_annealing_schedule<simdLen, Point>(newPin, P, ratio, C, parameters.frac, N, walk_length, error, a_vals, inv_covariance_matrix, rng);
+
+#ifdef VOLESTI_DEBUG
+    std::cout<<"All the variances of schedule_annealing computed in = "
+            << (double)clock()/(double)CLOCKS_PER_SEC-tstart2<<" sec"<<std::endl;
+    auto j=0;
+    for (auto avalIt = a_vals.begin(); avalIt!=a_vals.end(); avalIt++, j++){
+        std::cout<<"a_"<<j<<" = "<<*avalIt<<" ";
+    }
+    std::cout<<std::endl<<std::endl;
+#endif
+
+    // Initialization for the approximation of the ratios
+    unsigned int W = parameters.W;
+    unsigned int mm = a_vals.size()-1;
+    std::vector<NT> last_W2(W,0);
+    std::vector<NT> fn(mm,0);
+    std::vector<NT> its(mm,0);
+    VT lamdas;
+    lamdas.setZero(m);
+    
+    MT scaled_matrix = a_vals[0] * inv_covariance_matrix;
+    NT det_scaled = scaled_matrix.determinant();
+    NT vol = std::pow(M_PI, n / 2.0) / std::sqrt(det_scaled);
+
+    unsigned int i=0;
+
+    typedef typename std::vector<NT>::iterator viterator;
+    viterator itsIt = its.begin();
+    auto avalsIt = a_vals.begin();
+    viterator minmaxIt;
+
+
+#ifdef VOLESTI_DEBUG
+    std::cout<<"volume of the first gaussian = "<<vol<<"\n"<<std::endl;
+    std::cout<<"computing ratios..\n"<<std::endl;
+#endif
+
+    //iterate over the number of ratios
+    for (viterator fnIt = fn.begin();
+         fnIt != fn.end();
+         fnIt++, itsIt++, avalsIt++, i++)
+    {
+        //initialize convergence test
+        bool done = false;
+        NT curr_eps = error/std::sqrt((NT(mm)));
+        NT min_val = std::numeric_limits<NT>::min();
+        NT max_val = std::numeric_limits<NT>::max();
+        unsigned int min_index = W-1;
+        unsigned int max_index = W-1;
+        unsigned int index = 0;
+        unsigned int min_steps = 0;
+        std::vector<NT> last_W = last_W2;
+
+        // Set the radius for the ball walk
+        //creating the walk object
+        int dimension = newPin.dimension();
+        func_params f_params = func_params(Point(dimension), *avalsIt, -1, inv_covariance_matrix);
+
+        Func f(f_params);
+        Grad g(f_params);
+
+        ZeroFunctor<Point> zerof;
+
+        //create the crhmc problem
+        Input input = convert2crhmc_input<Input, Polytope, Func, Grad, ZeroFunctor<Point>>(newPin, f, g, zerof);
+
+        typedef crhmc_problem<Point, Input> CrhmcProblem;
+        CrhmcProblem problem = CrhmcProblem(input);
+        Point p = problem.center;
+
+        if(problem.terminate){return 0;}
+        problem.options.simdLen=simdLen;
+
+        //create the walk and do the burnIn
+        crhmc_walk_params params(input.df, p.dimension(), problem.options);
+        CRHMCWalkType walk = CRHMCWalkType(problem, p, input.df, input.f, params);
+
+        burnIn<CRHMCWalkType, crhmc_walk_params, simdLen, Grad, Func, CrhmcProblem>(
+            p, walk_length, rng, g, f, params, problem, walk, 2.5 * N);
+
+        while (!done || (*itsIt) < min_steps)
+        {
+            walk.apply(rng, walk_length);
+            p = walk.getPoint();
+            *itsIt = *itsIt + 1.0;
+            
+            *fnIt += eval_exp(p, inv_covariance_matrix, *(avalsIt+1), *avalsIt);
+
+            NT val = (*fnIt) / (*itsIt);
+
+            last_W[index] = val;
+            if (val <= min_val)
+            {
+                min_val = val;
+                min_index = index;
+            } else if (min_index == index)
+            {
+                minmaxIt = std::min_element(last_W.begin(), last_W.end());
+                min_val = *minmaxIt;
+                min_index = std::distance(last_W.begin(), minmaxIt);
+            }
+
+            if (val >= max_val)
+            {
+                max_val = val;
+                max_index = index;
+            } else if (max_index == index)
+            {
+                minmaxIt = std::max_element(last_W.begin(), last_W.end());
+                max_val = *minmaxIt;
+                max_index = std::distance(last_W.begin(), minmaxIt);
+            }
+
+            if ((max_val - min_val) / max_val <= curr_eps / 2.0)
+            {
+                done = true;
+            }
+
+            index = index % W + 1;
+            if (index == W) index = 0;
+        }
+#ifdef VOLESTI_DEBUG
+        std::cout << "ratio " << i << " = " << (*fnIt) / (*itsIt)
+                  << " N_" << i << " = " << *itsIt << std::endl;
+#endif
+        vol *= ((*fnIt) / (*itsIt));
+    }
+
+#ifdef VOLESTI_DEBUG
+        NT sum_of_steps = 0.0;
+        for(viterator it = its.begin(); it != its.end(); ++it) {
+            sum_of_steps += *it;
+        }
+        auto steps= int(sum_of_steps);
+        std::cout<<"\nTotal number of steps = "<<steps<<"\n"<<std::endl;
+#endif
+
+    return vol;
+}
+
+#endif // VOLUME_COOLING_NON_SPHERICAL_GAUSSIANS_CRHMC_HPP
\ No newline at end of file