diff --git a/examples/crhmc_cooling_gaussians/CMakeLists.txt b/examples/crhmc_cooling_gaussians/CMakeLists.txt
new file mode 100644
index 000000000..cd1a5a742
--- /dev/null
+++ b/examples/crhmc_cooling_gaussians/CMakeLists.txt
@@ -0,0 +1,132 @@
+# 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
+
+project( VolEsti )
+
+
+CMAKE_MINIMUM_REQUIRED(VERSION 3.11)
+
+set(CMAKE_ALLOW_LOOSE_LOOP_CONSTRUCTS true)
+
+# Locate Intel MKL root (in case it is enabled)
+if (APPLE)
+	set(MKLROOT /opt/intel/oneapi/mkl/latest)
+elseif(UNIX)
+  #set(MKLROOT /opt/intel/oneapi/mkl/latest)
+	set(MKLROOT $ENV{HOME}/intel/mkl)
+endif()
+
+if(COMMAND cmake_policy)
+  cmake_policy(SET CMP0003 NEW)
+endif(COMMAND cmake_policy)
+
+
+option(DISABLE_NLP_ORACLES "Disable non-linear oracles (used in collocation)" ON)
+option(BUILTIN_EIGEN "Use eigen from ../../external" OFF)
+option(USE_MKL "Use MKL library to build eigen" OFF)
+
+
+if(DISABLE_NLP_ORACLES)
+  add_definitions(-DDISABLE_NLP_ORACLES)
+else()
+  find_library(IFOPT NAMES libifopt_core.so PATHS /usr/local/lib)
+  find_library(IFOPT_IPOPT NAMES libifopt_ipopt.so PATHS /usr/local/lib)
+  find_library(GMP NAMES libgmp.so PATHS /usr/lib/x86_64-linux-gnu /usr/lib/i386-linux-gnu)
+  find_library(MPSOLVE NAMES libmps.so PATHS /usr/local/lib)
+  find_library(PTHREAD NAMES libpthread.so PATHS /usr/lib/x86_64-linux-gnu /usr/lib/i386-linux-gnu)
+  find_library(FFTW3 NAMES libfftw3.so.3 PATHS /usr/lib/x86_64-linux-gnu /usr/lib/i386-linux-gnu)
+
+  if (NOT IFOPT)
+
+    message(FATAL_ERROR "This program requires the ifopt library, and will not be compiled.")
+
+  elseif (NOT GMP)
+
+    message(FATAL_ERROR "This program requires the gmp library, and will not be compiled.")
+
+  elseif (NOT MPSOLVE)
+
+    message(FATAL_ERROR "This program requires the mpsolve library, and will not be compiled.")
+
+  elseif (NOT FFTW3)
+
+    message(FATAL_ERROR "This program requires the fftw3 library, and will not be compiled.")
+
+  else()
+    message(STATUS "Library ifopt found: ${IFOPT}")
+    message(STATUS "Library gmp found: ${GMP}")
+    message(STATUS "Library mpsolve found: ${MPSOLVE}")
+    message(STATUS "Library fftw3 found:" ${FFTW3})
+
+  endif(NOT IFOPT)
+
+endif(DISABLE_NLP_ORACLES)
+
+include("../../external/cmake-files/Eigen.cmake")
+GetEigen()
+
+include("../../external/cmake-files/Boost.cmake")
+GetBoost()
+
+include("../../external/cmake-files/LPSolve.cmake")
+GetLPSolve()
+
+include("../../external/cmake-files/QD.cmake")
+GetQD()
+
+# Find lpsolve library
+find_library(LP_SOLVE NAMES liblpsolve55.so PATHS /usr/lib/lp_solve)
+
+if (NOT LP_SOLVE)
+  message(FATAL_ERROR "This program requires the lp_solve library, and will not be compiled.")
+else ()
+  message(STATUS "Library lp_solve found: ${LP_SOLVE}")
+
+  set(CMAKE_EXPORT_COMPILE_COMMANDS "ON")
+
+  if (USE_MKL)
+    find_library(BLAS NAMES libblas.so libblas.dylib PATHS /usr/local/Cellar/lapack/3.9.1_1/lib /usr/lib/x86_64-linux-gnu /usr/lib/i386-linux-gnu /usr/local/Cellar/openblas/0.3.15_1/lib /usr/lib)
+    find_library(GFORTRAN NAME libgfortran.dylib PATHS /usr/local/Cellar/gcc/10.2.0_4/lib/gcc/10)
+    find_library(LAPACK NAME liblapack.dylib PATHS /usr/lib)
+    find_library(OPENMP NAME libiomp5.dylib PATHS /opt/intel/oneapi/compiler/2021.1.1/mac/compiler/lib)
+
+    include_directories (BEFORE ${MKLROOT}/include)
+    set(PROJECT_LIBS ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES} ${GFORTRAN_LIBRARIES})
+    set(MKL_LINK "-L${MKLROOT}/lib -Wl,-rpath,${MKLROOT}/lib -lmkl_intel_ilp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl")
+    add_definitions(-DEIGEN_USE_MKL_ALL)
+  else()
+    set(MKL_LINK "")
+  endif(USE_MKL)
+
+  include_directories (BEFORE ../../external)
+  include_directories (BEFORE ../../external/minimum_ellipsoid)
+  include_directories (BEFORE ../../include/)
+
+  # for Eigen
+  if (${CMAKE_VERSION} VERSION_LESS "3.12.0")
+    add_compile_options(-D "EIGEN_NO_DEBUG")
+  else ()
+    add_compile_definitions("EIGEN_NO_DEBUG")
+  endif ()
+
+
+  add_definitions(${CMAKE_CXX_FLAGS} "-std=c++17")  # enable C++17 standard
+  set(ADDITIONAL_FLAGS "-march=native -DSIMD_LEN=0 -DTIME_KEEPING")
+  add_definitions(${CMAKE_CXX_FLAGS} "-O3 -DTIME_KEEPING"  ${ADDITIONAL_FLAGS})  # optimization of the compiler
+  #add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lgsl")
+  add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lm")
+  add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-ldl")
+  add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-DBOOST_NO_AUTO_PTR")
+
+  add_executable(volume_example volume_example.cpp)
+  target_link_libraries(volume_example QD_LIB ${MKL_LINK} ${LP_SOLVE})
+
+endif()
diff --git a/examples/crhmc_cooling_gaussians/volume_example.cpp b/examples/crhmc_cooling_gaussians/volume_example.cpp
new file mode 100644
index 000000000..1a063e7d2
--- /dev/null
+++ b/examples/crhmc_cooling_gaussians/volume_example.cpp
@@ -0,0 +1,60 @@
+// 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
+
+#include "generators/known_polytope_generators.h"
+#include "random_walks/random_walks.hpp"
+#include "volume/volume_cooling_gaussians_crhmc.hpp"
+
+#include <iostream>
+#include <fstream>
+#include "misc/misc.h"
+
+typedef double NT;
+typedef Cartesian<NT> Kernel;
+typedef typename Kernel::Point Point;
+typedef BoostRandomNumberGenerator<boost::mt11213b, double> RandomNumberGenerator;
+typedef HPolytope<Point> HPOLYTOPE;
+
+void calculateAndVerifyVolume(HPOLYTOPE& polytope) {
+    int walk_len = 10;
+    NT e = 0.1;
+
+    RandomNumberGenerator rng(polytope.dimension());
+
+    NT volume = volume_cooling_gaussians<HPOLYTOPE, RandomNumberGenerator>(polytope, rng, e, walk_len);
+
+    std::cout << "Volume " << volume << std::endl;
+}
+
+int main() {
+    
+    HPOLYTOPE simplex = generate_simplex<HPOLYTOPE>(2, false);
+    std::cout << std::endl << "Simplex: " << std::endl;
+    simplex.print();
+    calculateAndVerifyVolume(simplex);
+    
+    HPOLYTOPE cube = generate_cube<HPOLYTOPE>(3, false);
+    std::cout << std::endl << "Cube: " << std::endl;
+    cube.print();
+    calculateAndVerifyVolume(cube);
+
+    HPOLYTOPE cross = generate_cross<HPOLYTOPE>(3, false);
+    std::cout << std::endl << "Cross: " << std::endl;
+    cross.print();
+    calculateAndVerifyVolume(cross);
+ 
+    HPOLYTOPE birkhoff = generate_birkhoff<HPOLYTOPE>(3);
+    std::cout << std::endl << "Birkhoff: " << std::endl;
+    birkhoff.print();
+    calculateAndVerifyVolume(birkhoff);
+    
+    return 0;
+}
diff --git a/include/volume/volume_cooling_gaussians_crhmc.hpp b/include/volume/volume_cooling_gaussians_crhmc.hpp
new file mode 100644
index 000000000..779f5b588
--- /dev/null
+++ b/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