diff --git a/examples/crhmc_sampling/simple_crhmc.cpp b/examples/crhmc_sampling/simple_crhmc.cpp
index 8f07065b0..86b236bca 100644
--- a/examples/crhmc_sampling/simple_crhmc.cpp
+++ b/examples/crhmc_sampling/simple_crhmc.cpp
@@ -25,7 +25,9 @@
 #include "random.hpp"
 #include <vector>
 #include "random_walks/random_walks.hpp"
+#include "preprocess/analytic_center_linear_ineq.h"
 #include "generators/known_polytope_generators.h"
+#include "generators/h_polytopes_generator.h"
 #include "helper_functions.hpp"
 
 template <int simdLen>
diff --git a/examples/sample_points/CMakeLists.txt b/examples/sample_points/CMakeLists.txt
new file mode 100644
index 000000000..aff504633
--- /dev/null
+++ b/examples/sample_points/CMakeLists.txt
@@ -0,0 +1,32 @@
+project( VolEsti-cpp-example )
+
+CMAKE_MINIMUM_REQUIRED(VERSION 3.11)
+
+add_definitions(-DDISABLE_NLP_ORACLES)
+
+include("../../external/cmake-files/Eigen.cmake")
+GetEigen()
+if (${CMAKE_VERSION} VERSION_LESS "3.12.0")
+  add_compile_options(-D "EIGEN_NO_DEBUG")
+else ()
+  add_compile_definitions("EIGEN_NO_DEBUG")
+endif ()
+
+include("../../external/cmake-files/Boost.cmake")
+GetBoost()
+
+include("../../external/cmake-files/LPSolve.cmake")
+GetLPSolve()
+
+include("../../external/cmake-files/QD.cmake")
+GetQD()
+
+set(CMAKE_EXPORT_COMPILE_COMMANDS "ON")
+
+add_definitions(${CMAKE_CXX_FLAGS} "-g")
+
+include_directories (BEFORE ../../external)
+include_directories (BEFORE ../../include)
+
+add_executable (example simple_example.cpp)
+target_link_libraries(example PUBLIC lp_solve QD_LIB)
\ No newline at end of file
diff --git a/examples/sample_points/simple_example.cpp b/examples/sample_points/simple_example.cpp
new file mode 100644
index 000000000..27b496cb4
--- /dev/null
+++ b/examples/sample_points/simple_example.cpp
@@ -0,0 +1,375 @@
+#include <chrono>
+#include "Eigen/Eigen"
+#include "cartesian_geom/cartesian_kernel.h"
+#include "convex_bodies/hpolytope.h"
+
+#include "diagnostics/diagnostics.hpp"
+//todo make headers automomous
+//#include "diagnostics/effective_sample_size.hpp"
+//#include "diagnostics/print_diagnostics.hpp"
+
+#include "generators/known_polytope_generators.h"
+
+#include "random_walks/random_walks.hpp"
+
+#include "sampling/sample_points.hpp"
+#include "sampling/sampling.hpp"
+
+#include "volume/volume_cooling_balls.hpp"
+
+typedef Cartesian<double> Kernel;
+typedef typename Kernel::Point Point;
+typedef BoostRandomNumberGenerator<boost::mt19937, double> RNGType;
+typedef HPolytope<Point> HPolytopeType;
+
+using NT = double;
+using MT = Eigen::Matrix<NT,Eigen::Dynamic,Eigen::Dynamic>;
+using VT = Eigen::Matrix<NT,Eigen::Dynamic,1>;
+
+using AcceleratedBilliardWalkType = AcceleratedBilliardWalk::template Walk<HPolytopeType, RNGType>;
+
+auto accelerated_billiard_walk(HPolytopeType& HP, RNGType& rng, unsigned int walk_len, unsigned int num_points)
+{
+    typedef RandomPointGenerator<AcceleratedBilliardWalkType> Generator;
+
+    MT Points(HP.dimension(), num_points);
+    Point q(HP.dimension()); // origin
+    AcceleratedBilliardWalkType walk(HP, q, rng);
+    walk.template get_starting_point(HP, q, q, 10, rng);
+    for (unsigned int i=0; i<num_points; ++i) {
+        walk.template apply(HP, q, 1, rng);
+        Points.col(i) = q.getCoefficients();
+    }
+
+    return Points;
+}
+
+template <typename PolytopeOrProblem, typename Walk, typename Distribution>
+void sample_points_eigen_matrix(PolytopeOrProblem const& HP, Point const& q, Walk const& walk,
+                                Distribution const& distr, RNGType rng, int walk_len, int rnum,
+                                int nburns)
+{
+    MT samples(q.dimension(), rnum);
+
+    auto start = std::chrono::steady_clock::now();
+
+    sample_points<64>(HP, q, walk, distr, rng, walk_len, rnum, nburns, samples);
+
+    auto end = std::chrono::steady_clock::now();
+    auto time = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
+    std::cout << " time= " << time;
+
+    // sample stats
+    unsigned int min_ess;
+    auto score = effective_sample_size<NT, VT, MT>(samples, min_ess);
+    std::cout << "\t ess= "  << min_ess << std::endl;
+    //print_diagnostics<NT, VT, MT>(samples, min_ess, std::cerr);
+}
+
+struct CustomFunctor {
+
+  // Custom density with neg log prob equal to c^T x
+  template
+  <
+      typename NT,
+      typename Point
+  >
+  struct parameters {
+    unsigned int order;
+    NT L; // Lipschitz constant for gradient
+    NT m; // Strong convexity constant
+    NT kappa; // Condition number
+    Point x0;
+
+    parameters(Point x0_) : order(2), L(1), m(1), kappa(1), x0(x0_) {};
+
+  };
+
+  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_) {};
+
+    // The index i represents the state vector index
+    Point operator() (unsigned int const& i, pts const& xs, NT const& t) const {
+      if (i == params.order - 1) {
+        Point y = (-1.0) * (xs[0] - params.x0);
+        return y;
+      } else {
+        return xs[i + 1]; // returns derivative
+      }
+    }
+  };
+
+  template<typename Point>
+  struct FunctionFunctor {
+    typedef typename Point::FT NT;
+
+    parameters<NT, Point> &params;
+
+    FunctionFunctor(parameters<NT, Point> &params_) : params(params_) {};
+
+    // The index i represents the state vector index
+    NT operator() (Point const& x) const {
+      Point y = x - params.x0;
+      return 0.5 * y.dot(y);
+    }
+  };
+};
+
+inline bool exists_check(const std::string &name) {
+  std::ifstream f(name.c_str());
+  return f.good();
+}
+
+template< typename SpMat, typename VT>
+void load_crhmc_problem(SpMat &A, VT &b, VT &lb, VT &ub, int &dimension,
+                        std::string problem_name) {
+   {
+    std::string fileName("./data/");
+    fileName.append(problem_name);
+    fileName.append(".mm");
+    if(!exists_check(fileName)){
+      std::cerr<<"Problem does not exist.\n";
+      exit(1);}
+    SpMat X;
+    loadMarket(X, fileName);
+    int m = X.rows();
+    dimension = X.cols() - 1;
+    A = X.leftCols(dimension);
+    b = VT(X.col(dimension));
+  }
+  {
+    std::string fileName("./data/");
+    fileName.append(problem_name);
+    fileName.append("_bounds.mm");
+    if(!exists_check(fileName)){
+      std::cerr<<"Problem does not exist.\n";
+      exit(1);}
+    SpMat bounds;
+    loadMarket(bounds, fileName);
+    lb = VT(bounds.col(0));
+    ub = VT(bounds.col(1));
+  }
+}
+
+
+int main() {
+
+    // 1. Inputs:
+
+    int d = 10;
+
+    // Generating a 3-dimensional cube centered at origin
+    HPolytopeType HP = generate_cube<HPolytopeType>(d, false);
+    std::cout<<"Polytope: \n";
+    HP.ComputeInnerBall();
+    //HP.print();
+    //std::cout<<"\n";
+
+    // Setup parameters for sampling
+    Point q(HP.dimension());
+    RNGType rng(HP.dimension());
+
+    // Generating a sparse polytope/problem
+    using SpMat = Eigen::SparseMatrix<NT>;
+    using ConstraintProblem =constraint_problem<SpMat, Point>;
+    std::string problem_name("simplex3");
+    std::cerr << "CRHMC on " << problem_name << "\n";
+    SpMat As;
+    VT b, lb, ub;
+    int dimension;
+    load_crhmc_problem(As, b, lb, ub, dimension, problem_name);
+    ConstraintProblem problem = ConstraintProblem(dimension);
+    problem.set_equality_constraints(As, b);
+    problem.set_bounds(lb, ub);
+
+
+    // 2. Walks
+
+    AcceleratedBilliardWalk abill_walk;
+    AcceleratedBilliardWalk abill_walk_custom(10); //user defined walk parameters
+    BallWalk ball_walk;
+    BilliardWalk bill_walk;
+    CDHRWalk cdhr_walk;
+    DikinWalk dikin_walk;
+    JohnWalk john_walk;
+    RDHRWalk rdhr_walk;
+    VaidyaWalk vaidya_walk;
+
+    GaussianBallWalk gball_walk;
+    GaussianCDHRWalk gcdhr_walk;
+    GaussianRDHRWalk grdhr_walk;
+    GaussianHamiltonianMonteCarloExactWalk ghmc_walk;
+
+    GaussianAcceleratedBilliardWalk gbill_walk;
+
+    ExponentialHamiltonianMonteCarloExactWalk ehmc_walk;
+
+    HamiltonianMonteCarloWalk hmc_walk;
+    NutsHamiltonianMonteCarloWalk nhmc_walk;
+    UnderdampedLangevinWalk uld_walk;
+    CRHMCWalk crhmc_walk;
+
+    // 3. Distributions
+
+    // Uniform
+    UniformDistribution udistr{};
+
+    // Spherical Gaussian
+    SphericalGaussianDistribution sgdistr{};
+
+    MT A(2, 2);
+    A << 0.25, 0.75,
+         0.75, 3.25;
+    Ellipsoid<Point> ell(A);    // origin centered ellipsoid
+    GaussianDistribution gdistr(ell);
+
+    // Exponential
+    NT variance = 1.0;
+    auto c = GetDirection<Point>::apply(HP.dimension(), rng, false);
+    ExponentialDistribution edistr(c, variance);
+
+    // LogConcave
+
+    std::pair<Point, NT> inner_ball = HP.ComputeInnerBall();
+    Point x0 = inner_ball.first;
+
+    // Reflective HMC and Remmannian HMC are using slightly different functor interfaces
+    // TODO: check if this could be unified
+
+    using NegativeGradientFunctorR = CustomFunctor::GradientFunctor<Point>;
+    using NegativeLogprobFunctorR = CustomFunctor::FunctionFunctor<Point>;
+    CustomFunctor::parameters<NT, Point> params_r(x0);
+    NegativeGradientFunctorR gr(params_r);
+    NegativeLogprobFunctorR fr(params_r);
+    LogConcaveDistribution logconcave_reflective(gr, fr, params_r.L);
+
+    using NegativeGradientFunctor = GaussianFunctor::GradientFunctor<Point>;
+    using NegativeLogprobFunctor = GaussianFunctor::FunctionFunctor<Point>;
+    using HessianFunctor = GaussianFunctor::HessianFunctor<Point>;
+    GaussianFunctor::parameters<NT, Point> params(x0, 0.5, 1);
+    NegativeGradientFunctor g(params);
+    NegativeLogprobFunctor f(params);
+    HessianFunctor h(params);
+    LogConcaveDistribution logconcave_crhmc(g, f, h, params.L);
+
+    LogConcaveDistribution logconcave_ref_gaus(g, f, params.L);
+
+    ZeroScalarFunctor<Point> f0;
+    ZeroFunctor<Point> g0;
+    ZeroFunctor<Point> h0;
+    LogConcaveDistribution logconcave_uniform(g0, f0, h0, params.L);
+    //LogConcaveDistribution logconcave_uniform2(g0, f0, params.L);
+
+
+
+    // 4. Sampling
+
+    using NT = double;
+    using MT = Eigen::Matrix<NT,Eigen::Dynamic,Eigen::Dynamic>;
+    using VT = Eigen::Matrix<NT,Eigen::Dynamic,1>;
+
+    int rnum = 100;
+    int nburns = 5;
+    int walk_len = 1;
+
+    MT samples(HP.dimension(), rnum);
+
+    // 1. the eigen matrix interface
+    std::cout << "uniform" << std::endl;
+    sample_points_eigen_matrix(HP, q, abill_walk, udistr, rng, 1, rnum, nburns);
+    sample_points_eigen_matrix(HP, q, abill_walk_custom, udistr, rng, walk_len, rnum, nburns);
+    sample_points_eigen_matrix(HP, q, ball_walk, udistr, rng, 2*d, rnum, nburns);
+    sample_points_eigen_matrix(HP, q, cdhr_walk, udistr, rng, 2*d, rnum, nburns);
+    sample_points_eigen_matrix(HP, q, dikin_walk, udistr, rng, walk_len, rnum, nburns);
+    sample_points_eigen_matrix(HP, q, john_walk, udistr, rng, walk_len, rnum, nburns);
+    sample_points_eigen_matrix(HP, q, rdhr_walk, udistr, rng, 2*d, rnum, nburns);
+    sample_points_eigen_matrix(HP, q, vaidya_walk, udistr, rng, walk_len, rnum, nburns);
+
+
+    std::cout << "shperical gaussian" << std::endl;
+    sample_points_eigen_matrix(HP, q, gball_walk, sgdistr, rng, walk_len, rnum, nburns);
+    sample_points_eigen_matrix(HP, q, gcdhr_walk, sgdistr, rng, walk_len, rnum, nburns);
+    sample_points_eigen_matrix(HP, q, grdhr_walk, sgdistr, rng, walk_len, rnum, nburns);
+    sample_points_eigen_matrix(HP, q, ghmc_walk, sgdistr, rng, walk_len, rnum, nburns);
+
+    std::cout << "general gaussian" << std::endl;
+    sample_points_eigen_matrix(HP, q, gbill_walk, gdistr, rng, walk_len, rnum, nburns);
+
+    std::cout << "exponential" << std::endl;
+    sample_points_eigen_matrix(HP, q, ehmc_walk, edistr, rng, walk_len, rnum, nburns);
+
+    std::cout << "logconcave" << std::endl;
+    sample_points_eigen_matrix(HP, q, hmc_walk, logconcave_reflective, rng, walk_len, rnum, nburns);
+    sample_points_eigen_matrix(HP, q, nhmc_walk, logconcave_reflective, rng, walk_len, rnum, nburns);
+    sample_points_eigen_matrix(HP, q, nhmc_walk, logconcave_ref_gaus, rng, walk_len, rnum, nburns);
+    sample_points_eigen_matrix(HP, q, uld_walk, logconcave_ref_gaus, rng, walk_len, rnum, nburns);
+
+    sample_points_eigen_matrix(HP, q, crhmc_walk, logconcave_crhmc, rng, walk_len, rnum, nburns);
+    // The following will compile but segfauls since walk and distribution are not compatible
+    //sample_points_eigen_matrix(HP, q, nhmc_walk, logconcave_crhmc, rng, walk_len, rnum, nburns);
+    sample_points_eigen_matrix(problem, q, crhmc_walk, logconcave_crhmc, rng, walk_len, rnum, nburns);
+    sample_points_eigen_matrix(problem, q, crhmc_walk, logconcave_uniform, rng, walk_len, 100000, nburns);
+    sample_points_eigen_matrix(HP, q, crhmc_walk, logconcave_uniform, rng, walk_len, rnum, nburns);
+    //sample_points_eigen_matrix(problem, q, nhmc_walk, logconcave_uniform2, rng, walk_len, rnum, nburns);
+
+    using NT_fromMT = VT::Scalar;
+    NT_fromMT test;
+
+    std::cout << "fix the following" << std::endl;
+    // TODO: fix
+    // Does not converge because of the starting point
+    // Also ess returns rnum instead of 0
+    //sample_points_eigen_matrix(HP, q, bill_walk, udistr, rng, walk_len, rnum, nburns);
+
+    // Does not compile because of walk-distribution combination
+    //sample_points_eigen_matrix(HP, q, abill_walk, gdistr, rng, walk_len, rnum, nburns);
+/*
+    std::cout << "std::vector interface" << std::endl;
+    // 2. the std::vector interface
+    std::vector<Point> points;
+    sample_points(HP, q, cdhr_walk, udistr, rng, walk_len, rnum, nburns, points);
+    for (auto& point : points)
+    {
+    //    std::cout << point.getCoefficients().transpose() << "\n";
+    }
+*/
+/*
+    // 3. the old interface
+    // different billiard walks
+    typedef BilliardWalk::template Walk<HPolytopeType, RNGType> BilliardWalkType;
+    typedef AcceleratedBilliardWalk::template Walk<HPolytopeType, RNGType> AcceleratedBilliardWalkType;
+    typedef RandomPointGenerator<AcceleratedBilliardWalkType> Generator;
+    std::vector<Point> randPoints;
+    PushBackWalkPolicy push_back_policy;
+    Generator::apply(HP, q, rnum, walk_len, randPoints, push_back_policy, rng);
+    for (auto& point : randPoints)
+    {
+    //    std::cout << point.getCoefficients().transpose() << "\n";
+    }
+*/
+/*
+    unsigned int walkL = 10, numpoints = 10000, nburns = 0, d = HP.dimension();
+    Point StartingPoint(d);
+    std::list<Point> randPoints;
+
+//    gaussian_sampling<AcceleratedBilliardWalk>(points, HP, rng, walkL, numpoints, 1.0,
+    //                               StartingPoint, nburns);
+
+    double variance = 1.0;
+
+    Point c(HP.dimension());
+    HP.set_InnerBall(std::pair<Point,double>(Point(HP.dimension()), 1.0));
+    c = GetDirection<Point>::apply(HP.dimension(), rng, false);
+    //ExponentialHamiltonianMonteCarloExactWalk
+    exponential_sampling<NutsHamiltonianMonteCarloWalk>(points, HP, rng, walkL, numpoints, c, variance,
+                                StartingPoint, nburns);
+
+*/
+    return 0;
+}
\ No newline at end of file
diff --git a/include/preprocess/inscribed_ellipsoid_rounding.hpp b/include/preprocess/inscribed_ellipsoid_rounding.hpp
index 4b027f61b..cd99a519c 100644
--- a/include/preprocess/inscribed_ellipsoid_rounding.hpp
+++ b/include/preprocess/inscribed_ellipsoid_rounding.hpp
@@ -37,7 +37,7 @@ compute_inscribed_ellipsoid(Custom_MT A, VT b, VT const& x0,
     return {};
 }
 
-template 
+template
 <
     typename MT,
     typename VT,
@@ -45,7 +45,7 @@ template
     typename Polytope,
     int ellipsoid_type = EllipsoidType::MAX_ELLIPSOID
 >
-std::tuple<MT, VT, NT> inscribed_ellipsoid_rounding(Polytope &P, 
+std::tuple<MT, VT, NT> inscribed_ellipsoid_rounding(Polytope &P,
                                                     unsigned int const max_iterations = 5,
                                                     NT const max_eig_ratio = NT(6))
 {
@@ -54,7 +54,7 @@ std::tuple<MT, VT, NT> inscribed_ellipsoid_rounding(Polytope &P,
     return inscribed_ellipsoid_rounding<MT, VT, NT>(P, Point(x), max_iterations, max_eig_ratio);
 }
 
-template 
+template
 <
     typename MT,
     typename VT,
@@ -63,7 +63,7 @@ template
     typename Point,
     int ellipsoid_type = EllipsoidType::MAX_ELLIPSOID
 >
-std::tuple<MT, VT, NT> inscribed_ellipsoid_rounding(Polytope &P, 
+std::tuple<MT, VT, NT> inscribed_ellipsoid_rounding(Polytope &P,
                                                     Point const& InnerPoint,
                                                     unsigned int const max_iterations = 5,
                                                     NT const max_eig_ratio = NT(6))
@@ -77,9 +77,9 @@ std::tuple<MT, VT, NT> inscribed_ellipsoid_rounding(Polytope &P,
     while (true)
     {
         // Compute the desired inscribed ellipsoid in P
-        std::tie(E, center, converged) = 
+        std::tie(E, center, converged) =
             compute_inscribed_ellipsoid<MT, ellipsoid_type>(P.get_mat(), P.get_vec(), x0, maxiter, tol, reg);
-        
+
         E = (E + E.transpose()) / 2.0;
         E += MT::Identity(d, d)*std::pow(10, -8.0); //normalize E
 
@@ -89,7 +89,7 @@ std::tuple<MT, VT, NT> inscribed_ellipsoid_rounding(Polytope &P,
         // Computing eigenvalues of E
         Spectra::DenseSymMatProd<NT> op(E);
         // The value of ncv is chosen empirically
-        Spectra::SymEigsSolver<NT, Spectra::SELECT_EIGENVALUE::BOTH_ENDS, 
+        Spectra::SymEigsSolver<NT, Spectra::SELECT_EIGENVALUE::BOTH_ENDS,
                                Spectra::DenseSymMatProd<NT>> eigs(&op, 2, std::min(std::max(10, int(d)/5), int(d)));
         eigs.init();
         int nconv = eigs.compute();
@@ -128,4 +128,4 @@ std::tuple<MT, VT, NT> inscribed_ellipsoid_rounding(Polytope &P,
     return result;
 }
 
-#endif 
+#endif
diff --git a/include/random_walks/compute_diameter.hpp b/include/random_walks/compute_diameter.hpp
index 71007d9e4..641cae6dd 100644
--- a/include/random_walks/compute_diameter.hpp
+++ b/include/random_walks/compute_diameter.hpp
@@ -24,10 +24,7 @@
 
 template <typename GenericPolytope>
 struct compute_diameter
-{
-    template <typename NT>
-    static NT compute(GenericPolytope) {return NT(0);}
-};
+{};
 
 
 template <typename Point, typename MT>
@@ -40,6 +37,16 @@ struct compute_diameter<HPolytope<Point, MT>>
     }
 };
 
+template <typename Point>
+struct compute_diameter<const HPolytope<Point>>
+{
+    template <typename NT>
+    static NT compute(HPolytope<Point> const& P)
+    {
+        return NT(2) * std::sqrt(NT(P.dimension())) * P.InnerBall().second;
+    }
+};
+
 template <typename Point>
 struct compute_diameter<Spectrahedron<Point>>
 {
diff --git a/include/random_walks/crhmc/crhmc_walk.hpp b/include/random_walks/crhmc/crhmc_walk.hpp
index 9f1dd24af..3bffa4553 100644
--- a/include/random_walks/crhmc/crhmc_walk.hpp
+++ b/include/random_walks/crhmc/crhmc_walk.hpp
@@ -43,6 +43,16 @@ struct CRHMCWalk {
       eta = 1.0 / (dim * sqrt(F.params.L));
       momentum = 1 - std::min(1.0, eta / effectiveStepSize);
     }
+    parameters(NT const& L,
+               unsigned int dim,
+               Opts &user_options,
+               NT epsilon_ = 2)
+      : options(user_options)
+    {
+      epsilon = epsilon_;
+      eta = 1.0 / (dim * sqrt(L));
+      momentum = 1 - std::min(1.0, eta / effectiveStepSize);
+    }
   };
 
   template
diff --git a/include/random_walks/gaussian_cdhr_walk.hpp b/include/random_walks/gaussian_cdhr_walk.hpp
index 04ebedbc4..e267f368c 100644
--- a/include/random_walks/gaussian_cdhr_walk.hpp
+++ b/include/random_walks/gaussian_cdhr_walk.hpp
@@ -87,7 +87,7 @@ struct Walk
          Point const& p,
          NT const& a_i,
          RandomNumberGenerator& rng,
-         parameters&)
+         parameters const&)
     {
         initialize(P, p, a_i, rng);
     }
diff --git a/include/random_walks/gaussian_rdhr_walk.hpp b/include/random_walks/gaussian_rdhr_walk.hpp
index b7a53abe7..29cb8ea7f 100644
--- a/include/random_walks/gaussian_rdhr_walk.hpp
+++ b/include/random_walks/gaussian_rdhr_walk.hpp
@@ -84,7 +84,7 @@ struct Walk
     {}
 
     Walk(Polytope&, Point const&, NT const&, RandomNumberGenerator&,
-         parameters&)
+         parameters const&)
     {}
 
     template
diff --git a/include/random_walks/hamiltonian_monte_carlo_walk.hpp b/include/random_walks/hamiltonian_monte_carlo_walk.hpp
index 23342ea57..f696c4b1a 100644
--- a/include/random_walks/hamiltonian_monte_carlo_walk.hpp
+++ b/include/random_walks/hamiltonian_monte_carlo_walk.hpp
@@ -40,6 +40,12 @@ struct HamiltonianMonteCarloWalk {
       // eta = 1.0 /
       //   (sqrt(20 * F.params.L * pow(dim, 3)));
     }
+
+    parameters(NT const& L, unsigned int dim, NT epsilon_ = 2)
+    {
+      epsilon = epsilon_;
+      eta = 1.0 / (dim * sqrt(L));
+    }
   };
 
   template
diff --git a/include/random_walks/langevin_walk.hpp b/include/random_walks/langevin_walk.hpp
index df23dbacf..1e2a9c786 100644
--- a/include/random_walks/langevin_walk.hpp
+++ b/include/random_walks/langevin_walk.hpp
@@ -44,6 +44,13 @@ struct UnderdampedLangevinWalk {
       //                       pow(epsilon, 2.0 / 3) *
       //                       pow(log(1.0 / epsilon), - 1.0 / 3));
     }
+
+    parameters(NT const& L, unsigned int dim, NT epsilon_ = 1e-4)
+    {
+      epsilon = epsilon_;
+      u = 1.0 / L;
+      eta = 1.0 / (sqrt(20 * L));
+    }
   };
 
   template
diff --git a/include/random_walks/nuts_hmc_walk.hpp b/include/random_walks/nuts_hmc_walk.hpp
index 857064c0e..e6408864f 100644
--- a/include/random_walks/nuts_hmc_walk.hpp
+++ b/include/random_walks/nuts_hmc_walk.hpp
@@ -8,7 +8,7 @@
 // Licensed under GNU LGPL.3, see LICENCE file
 
 // References
-// Matthew D. Hoffman, Andrew Gelman. "The No-U-Turn Sampler: 
+// Matthew D. Hoffman, Andrew Gelman. "The No-U-Turn Sampler:
 // Adaptively Setting Path Lengths in Hamiltonian Monte Carlo", 2011.
 
 // Comment: Compared to [Matthew D. Hoffman, Andrew Gelman, 2011]
@@ -42,6 +42,12 @@ struct NutsHamiltonianMonteCarloWalk {
       epsilon = epsilon_;
       eta = F.params.L > 0 ? 10.0 / (dim * sqrt(F.params.L)) : 0.005;
     }
+
+    parameters(NT const& L, unsigned int dim, NT epsilon_ = 2)
+    {
+      epsilon = epsilon_;
+      eta = L > 0 ? 10.0 / (dim * sqrt(L)) : 0.005;
+    }
   };
 
   template
@@ -149,7 +155,7 @@ struct NutsHamiltonianMonteCarloWalk {
       Point p = x;
       NT L;
 
-      if ((solver->get_bounds())[0] == NULL) 
+      if ((solver->get_bounds())[0] == NULL)
       {
         L = (NT(100) / NT(dim)) * (NT(100) / NT(dim));
       }
@@ -181,12 +187,12 @@ struct NutsHamiltonianMonteCarloWalk {
 
       // Pick a random velocity
       v = GetDirection<Point>::apply(dim, rng, false);
-      
+
       v_pl = v;
       v_min = NT(-1) * v;
       X_pl = x;
       X_min = x;
-      
+
       NT h1 = hamiltonian(x, v);
 
       NT uu = std::log(rng.sample_urdist()) - h1;
@@ -203,14 +209,14 @@ struct NutsHamiltonianMonteCarloWalk {
       while (s)
       {
         j++;
-        
+
         if (burnin)
         {
           na = std::pow(NT(2), NT(j));
         }
 
         NT dir = rng.sample_urdist();
-                
+
         if (dir > 0.5)
         {
           v = v_pl;
@@ -222,7 +228,7 @@ struct NutsHamiltonianMonteCarloWalk {
           X = X_min;
         }
         X_rnd_j = X;
-                
+
         int x_counting = 0;
         int num_samples = int(std::pow(NT(2), NT(j)));
         accepted = false;
@@ -256,76 +262,76 @@ struct NutsHamiltonianMonteCarloWalk {
           }
 
           bool pos_state = false;
-          if (uu < -hj) 
+          if (uu < -hj)
           {
             pos_state = true;
             pos_state_single = true;
             x_counting = x_counting + 1;
             x_counting_total = x_counting_total + 1;
           }
-                 
-          if (k == 1) 
+
+          if (k == 1)
           {
-            if (dir > 0.5) 
+            if (dir > 0.5)
             {
               X_min_j = X;
               v_min_j = v;
-            } 
-            else 
+            }
+            else
             {
               X_pl_j = X;
               v_pl_j = v;
             }
           }
-          if (k == num_samples) 
+          if (k == num_samples)
           {
-            if (dir > 0.5) 
+            if (dir > 0.5)
             {
               x_pl_min = X - X_min_j;
-              if ((x_pl_min.dot(v) < 0) || (x_pl_min.dot(v_min_j) < 0)) 
+              if ((x_pl_min.dot(v) < 0) || (x_pl_min.dot(v_min_j) < 0))
               {
                 s = false;
               }
-            } 
-            else 
+            }
+            else
             {
               x_pl_min = X_pl_j - X;
-              if ((x_pl_min.dot(v) < 0) || (x_pl_min.dot(v_pl_j) < 0)) 
+              if ((x_pl_min.dot(v) < 0) || (x_pl_min.dot(v_pl_j) < 0))
               {
                 s = false;
               }
             }
           }
-          if ((rng.sample_urdist() < (1/NT(x_counting))) && pos_state) 
+          if ((rng.sample_urdist() < (1/NT(x_counting))) && pos_state)
           {
             X_rnd_j = X;
           }
         }
 
-        if (dir > 0.5) 
+        if (dir > 0.5)
         {
           X_pl = X;
           v_pl = v;
-        } 
-        else 
+        }
+        else
         {
           X_min = X;
           v_min = v;
         }
-                
-        if (s && (rng.sample_urdist() < (NT(x_counting) / NT(x_counting_total)))) 
+
+        if (s && (rng.sample_urdist() < (NT(x_counting) / NT(x_counting_total))))
         {
           x = X_rnd_j;
           if (pos_state_single)
-          { 
+          {
             updated = true;
           }
         }
-                
-        if (s) 
+
+        if (s)
         {
           x_pl_min = X_pl - X_min;
-          if ((x_pl_min.dot(v_min) < 0) || (x_pl_min.dot(v_pl) < 0)) 
+          if ((x_pl_min.dot(v_min) < 0) || (x_pl_min.dot(v_pl) < 0))
           {
               s = false;
           }
diff --git a/include/random_walks/uniform_accelerated_billiard_walk.hpp b/include/random_walks/uniform_accelerated_billiard_walk.hpp
index f91abee4b..daacda82b 100644
--- a/include/random_walks/uniform_accelerated_billiard_walk.hpp
+++ b/include/random_walks/uniform_accelerated_billiard_walk.hpp
@@ -181,7 +181,7 @@ struct AcceleratedBilliardWalk
         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 >; 
+        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>
@@ -250,7 +250,7 @@ struct AcceleratedBilliardWalk
 
                 it = 0;
                 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;
@@ -275,7 +275,7 @@ struct AcceleratedBilliardWalk
                 it++;
 
                 while (it < _rho)
-                {   
+                {
                     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,
@@ -334,7 +334,7 @@ struct AcceleratedBilliardWalk
         <
             typename GenericPolytope
         >
-        inline void parameters_burnin(GenericPolytope &P, 
+        inline void parameters_burnin(GenericPolytope &P,
                                      Point const& center,
                                      unsigned int const& num_points,
                                      unsigned int const& walk_length,
@@ -346,7 +346,7 @@ struct AcceleratedBilliardWalk
             pointset.push_back(_p);
             NT rad = NT(0), max_dist, Lmax = get_delta(), radius = P.InnerBall().second;
 
-            for (int i = 0; i < num_points; i++) 
+            for (int i = 0; i < num_points; i++)
             {
                 Point p = GetPointInDsphere<Point>::apply(P.dimension(), radius, rng);
                 p += center;
@@ -374,7 +374,7 @@ struct AcceleratedBilliardWalk
             pointset.clear();
         }
 
-        
+
 
         inline void update_delta(NT L)
         {
@@ -441,11 +441,11 @@ struct AcceleratedBilliardWalk
             }
         }
 
-        inline double get_max_distance(std::vector<Point> &pointset, Point const& q, double &rad) 
+        inline double get_max_distance(std::vector<Point> &pointset, Point const& q, double &rad)
         {
             double dis = -1.0, temp_dis;
             int jj = 0;
-            for (auto vecit = pointset.begin(); vecit!=pointset.end(); vecit++, jj++) 
+            for (auto vecit = pointset.begin(); vecit!=pointset.end(); vecit++, jj++)
             {
                 temp_dis = (q.getCoefficients() - (*vecit).getCoefficients()).norm();
                 if (temp_dis > dis) {
diff --git a/include/random_walks/uniform_dikin_walk.hpp b/include/random_walks/uniform_dikin_walk.hpp
index e791b2544..326413d2c 100644
--- a/include/random_walks/uniform_dikin_walk.hpp
+++ b/include/random_walks/uniform_dikin_walk.hpp
@@ -57,20 +57,22 @@ struct DikinWalk
         typedef typename Polytope::VT VT;
         typedef typename Point::FT NT;
 
-        Walk(Polytope &P, Point &p, RandomNumberGenerator &)
+        template <typename GenericPolytope>
+        Walk(GenericPolytope& P, Point const& p, RandomNumberGenerator &)
         {
             MT A = P.get_mat();
             VT b = P.get_vec(), _vec_point = VT::Zero(P.dimension()), p0 = p.getCoefficients();
-            NT r = P.ComputeInnerBall().second;
+            NT r = P.InnerBall().second;
             dikinw = DikinWalker<NT>(p0, A, b, r);
         }
 
-        Walk(Polytope &P, Point & p, RandomNumberGenerator &, parameters const& params)
+        template <typename GenericPolytope>
+        Walk(GenericPolytope& P, Point const& p, RandomNumberGenerator &, parameters const& params)
         {
             MT A = P.get_mat();
             VT b = P.get_vec(), _vec_point = VT::Zero(P.dimension()), p0 = p.getCoefficients();
             NT r = params.set_L ? params.m_L
-                          : P.ComputeInnerBall().second;
+                          : P.InnerBall().second;
             dikinw = DikinWalker<NT>(p0, A, b, r);
         }
 
diff --git a/include/random_walks/uniform_john_walk.hpp b/include/random_walks/uniform_john_walk.hpp
index 2c6bc42e3..cde92587a 100644
--- a/include/random_walks/uniform_john_walk.hpp
+++ b/include/random_walks/uniform_john_walk.hpp
@@ -49,20 +49,22 @@ struct JohnWalk
         typedef typename Polytope::VT VT;
         typedef typename Point::FT NT;
 
-        Walk(Polytope &P, Point &p, RandomNumberGenerator &)
+        template <typename GenericPolytope>
+        Walk(GenericPolytope& P, Point const& p, RandomNumberGenerator &)
         {
             MT A = P.get_mat();
             VT b = P.get_vec(), _vec_point = VT::Zero(P.dimension()), p0 = p.getCoefficients();
-            NT r = P.ComputeInnerBall().second;
+            NT r = P.InnerBall().second;
             johnw = JohnWalker<NT>(p0, A, b, r);
         }
 
-        Walk(Polytope &P, Point & p, RandomNumberGenerator &, parameters const& params)
+        template <typename GenericPolytope>
+        Walk(GenericPolytope& P, Point const& p, RandomNumberGenerator &, parameters const& params)
         {
             MT A = P.get_mat();
             VT b = P.get_vec(), _vec_point = VT::Zero(P.dimension()), p0 = p.getCoefficients();
             NT r = params.set_L ? params.m_L
-                          : P.ComputeInnerBall().second;
+                          : P.InnerBall().second;
             johnw = JohnWalker<NT>(p0, A, b, r);
         }
 
diff --git a/include/random_walks/uniform_vaidya_walk.hpp b/include/random_walks/uniform_vaidya_walk.hpp
index 5bbbe9912..fd076cc28 100644
--- a/include/random_walks/uniform_vaidya_walk.hpp
+++ b/include/random_walks/uniform_vaidya_walk.hpp
@@ -50,20 +50,22 @@ struct VaidyaWalk
         typedef typename Polytope::VT VT;
         typedef typename Point::FT NT;
 
-        Walk(Polytope &P, Point &p, RandomNumberGenerator &)
+        template <typename GenericPolytope>
+        Walk(GenericPolytope& P, Point const& p, RandomNumberGenerator &)
         {
             MT A = P.get_mat();
             VT b = P.get_vec(), _vec_point = VT::Zero(P.dimension()), p0 = p.getCoefficients();
-            NT r = P.ComputeInnerBall().second;
+            NT r = P.InnerBall().second;
             vaidyaw = VaidyaWalker<NT>(p0, A, b, r);
         }
 
-        Walk(Polytope &P, Point & p, RandomNumberGenerator &, parameters const& params)
+        template <typename GenericPolytope>
+        Walk(GenericPolytope& P, Point const& p, RandomNumberGenerator &, parameters const& params)
         {
             MT A = P.get_mat();
             VT b = P.get_vec(), _vec_point = VT::Zero(P.dimension()), p0 = p.getCoefficients();
             NT r = params.set_L ? params.m_L
-                          : P.ComputeInnerBall().second;
+                          : P.InnerBall().second;
             vaidyaw = VaidyaWalker<NT>(p0, A, b, r);
         }
 
diff --git a/include/sampling/sample_points.hpp b/include/sampling/sample_points.hpp
new file mode 100644
index 000000000..bc6d341cf
--- /dev/null
+++ b/include/sampling/sample_points.hpp
@@ -0,0 +1,347 @@
+// VolEsti (volume computation and sampling library)
+
+// Copyright (c) 2024 Vissarion Fisikopoulos
+
+// Licensed under GNU LGPL.3, see LICENCE file
+
+#ifndef SAMPLE_POINTS_HPP
+#define SAMPLE_POINTS_HPP
+
+#include "preprocess/crhmc/crhmc_input.h"
+#include "preprocess/crhmc/crhmc_problem.h"
+#include "sampling/sampling.hpp"
+
+template <typename NT = double>
+struct UniformDistribution
+{
+    using CartesianPoint = point<Cartesian<NT>>;
+    using Func = ZeroScalarFunctor<CartesianPoint>;
+    using Grad = ZeroFunctor<CartesianPoint>;
+    using Hess = ZeroFunctor<CartesianPoint>;
+};
+
+template <typename NT = double>
+struct SphericalGaussianDistribution
+{
+    SphericalGaussianDistribution()
+        :   variance(1.0)
+    {}
+    SphericalGaussianDistribution(NT _variance)
+        :   variance(_variance)
+    {}
+    NT variance;
+};
+
+template <typename NT = double>
+struct GaussianDistribution
+{
+    using CartesianEllipsoid = Ellipsoid<point<Cartesian<NT>>>;
+
+    GaussianDistribution() {}
+
+    GaussianDistribution(CartesianEllipsoid _ellipsoid)
+        :   ellipsoid(_ellipsoid)
+    {}
+    CartesianEllipsoid ellipsoid;
+};
+
+template <typename NT = double>
+struct ExponentialDistribution
+{
+    using CartesianPoint = point<Cartesian<NT>>;
+
+    ExponentialDistribution() {}
+
+    ExponentialDistribution(CartesianPoint _c, NT _T)
+        :   c(_c)
+        ,   T(_T)
+    {}
+    CartesianPoint c;
+    NT T;
+};
+
+template <typename NT = double>
+struct LogConcaveDistribution
+{
+    using CartesianPoint = point<Cartesian<NT>>;
+
+    template <typename GradientFunctor, typename FunctionFunctor>
+    LogConcaveDistribution(GradientFunctor g, FunctionFunctor f, NT _L)
+        : grad(g)
+        , func(f)
+        , L(_L)
+    {}
+
+    template <typename GradientFunctor, typename FunctionFunctor, typename HessianFunctor>
+    LogConcaveDistribution(GradientFunctor g, FunctionFunctor f, HessianFunctor h, NT _L)
+        : grad_point(g)
+        , func(f)
+        , hess(h)
+        , L(_L)
+    {}
+
+    std::function<CartesianPoint(unsigned int const&, std::vector<CartesianPoint> const&, NT const&)> grad;
+    std::function<CartesianPoint(CartesianPoint const&)> grad_point;
+    std::function<NT(CartesianPoint const&)> func;
+    std::function<CartesianPoint(CartesianPoint const&)> hess;
+    NT L;
+};
+
+namespace detail
+{
+
+template <typename T1, typename T2>
+struct AlwaysFalse : std::false_type {};
+
+template <typename WalkType, typename Distribution>
+struct sample_points
+{
+    static_assert(AlwaysFalse<WalkType, Distribution>::value, "Not implemented for this combination of walk and distribution.");
+};
+
+}
+
+template
+<
+    int simdLen = 1,
+    typename Polytope,
+    typename Point,
+    typename WalkType,
+    typename Distribution,
+    typename RandomNumberGenerator,
+    typename NT
+>
+void sample_points(Polytope& P, // TODO: make it a const&
+                   Point const& starting_point,
+                   WalkType const& walk_with_parameters,
+                   Distribution const& distribution,
+                   RandomNumberGenerator& rng,
+                   unsigned int const& walk_len,
+                   unsigned int const& rnum,
+                   unsigned int const& nburns,
+                   Eigen::Matrix<NT,Eigen::Dynamic,Eigen::Dynamic>& samples)
+{
+    if constexpr ((std::is_same<WalkType, AcceleratedBilliardWalk>::value
+                || std::is_same<WalkType, BallWalk>::value
+                || std::is_same<WalkType, BilliardWalk>::value
+                || std::is_same<WalkType, CDHRWalk>::value
+                || std::is_same<WalkType, DikinWalk>::value
+                || std::is_same<WalkType, JohnWalk>::value
+                || std::is_same<WalkType, RDHRWalk>::value
+                || std::is_same<WalkType, VaidyaWalk>::value)
+            && std::is_same<Distribution, UniformDistribution<NT>>::value)
+    {
+        typename WalkType::template Walk<Polytope, RandomNumberGenerator>
+            walk(P, starting_point, rng, walk_with_parameters.param);
+
+        Point p = starting_point;
+
+        for (unsigned int i = 0; i < nburns; ++i)
+        {
+            walk.apply(P, p, walk_len, rng);
+        }
+
+        for (unsigned int i = 0; i < rnum; ++i)
+        {
+            walk.apply(P, p, walk_len, rng);
+            samples.col(i) = p.getCoefficients();
+        }
+    }
+    else if constexpr ((std::is_same<WalkType, GaussianBallWalk>::value
+                     || std::is_same<WalkType, GaussianCDHRWalk>::value
+                     || std::is_same<WalkType, GaussianHamiltonianMonteCarloExactWalk>::value
+                     || std::is_same<WalkType, GaussianRDHRWalk>::value)
+                    && std::is_same<Distribution, SphericalGaussianDistribution<NT>>::value)
+    {
+        typename WalkType::template Walk<Polytope, RandomNumberGenerator>
+            walk(P, starting_point, distribution.variance, rng, walk_with_parameters.param);
+
+        Point p = starting_point;
+
+        for (unsigned int i = 0; i < nburns; ++i)
+        {
+            walk.apply(P, p, distribution.variance, walk_len, rng);
+        }
+
+        for (unsigned int i = 0; i < rnum; ++i)
+        {
+            walk.apply(P, p, distribution.variance, walk_len, rng);
+            samples.col(i) = p.getCoefficients();
+        }
+    }
+    else if constexpr (std::is_same<WalkType, GaussianAcceleratedBilliardWalk>::value
+                    && std::is_same<Distribution, GaussianDistribution<NT>>::value)
+    {
+        typename WalkType::template Walk<Polytope, RandomNumberGenerator>
+            walk(P, starting_point, distribution.ellipsoid, rng, walk_with_parameters.param);
+
+        Point p = starting_point;
+
+        for (unsigned int i = 0; i < nburns; ++i)
+        {
+            walk.apply(P, p, distribution.ellipsoid, walk_len, rng);
+        }
+
+        for (unsigned int i = 0; i < rnum; ++i)
+        {
+            walk.apply(P, p, distribution.ellipsoid, walk_len, rng);
+            samples.col(i) = p.getCoefficients();
+        }
+    }
+    else if constexpr (std::is_same<WalkType, ExponentialHamiltonianMonteCarloExactWalk>::value
+                    && std::is_same<Distribution, ExponentialDistribution<NT>>::value)
+    {
+        typename WalkType::template Walk<Polytope, RandomNumberGenerator>
+            walk(P, starting_point, distribution.c, distribution.T, rng, walk_with_parameters.param);
+
+        Point p = starting_point;
+
+        for (unsigned int i = 0; i < nburns; ++i)
+        {
+            walk.apply(P, p, walk_len, rng);
+        }
+
+        for (unsigned int i = 0; i < rnum; ++i)
+        {
+            walk.apply(P, p, walk_len, rng);
+            samples.col(i) = p.getCoefficients();
+        }
+    }
+    else if constexpr ((std::is_same<WalkType, HamiltonianMonteCarloWalk>::value
+                     || std::is_same<WalkType, NutsHamiltonianMonteCarloWalk>::value
+                     || std::is_same<WalkType, UnderdampedLangevinWalk>::value)
+                    && std::is_same<Distribution, LogConcaveDistribution<NT>>::value)
+    {
+        using HPolytope = typename std::remove_const<Polytope>::type;
+
+        using Solver = LeapfrogODESolver<Point, NT, HPolytope, decltype(distribution.grad)>;
+
+        std::vector<Point> xs;
+        unsigned int i = 0;
+        NT t = 1.0;
+
+        typename WalkType::template parameters
+        <
+            NT,
+            decltype(distribution.grad)
+        > hmc_params(distribution.L, P.dimension());
+
+        Point walk_p = starting_point; //TODO: avoid the copy
+        auto g = distribution.grad;
+        auto f = distribution.func;
+
+        HPolytope P_copy = P; //TODO: avoid the copy
+        typename WalkType::template Walk
+        <
+            Point, HPolytope, RandomNumberGenerator,
+            decltype(distribution.grad), decltype(distribution.func), Solver
+        >
+        walk(&P_copy, walk_p, g, f, hmc_params);
+
+        Point p = starting_point;
+
+        //TODO: burnin for nuts
+
+        for (int i = 0; i < rnum; i++)
+        {
+            walk.apply(rng, walk_len);
+            samples.col(i) = walk.x.getCoefficients();
+        }
+    }
+    else if constexpr ((std::is_same<WalkType, CRHMCWalk>::value)
+                    && std::is_same<Distribution, LogConcaveDistribution<NT>>::value)
+    {
+        using HPolytope = typename std::remove_const<Polytope>::type;
+        HPolytope HP = P; //TODO: avoid the copy
+
+        int dimension = HP.dimension();
+
+        auto G = distribution.grad_point;
+        auto F = distribution.func;
+        auto H = distribution.hess;
+
+        using NegativeLogprobFunctor = decltype(distribution.func);
+        using NegativeGradientFunctor = decltype(distribution.grad_point);
+        using HessianFunctor = decltype(distribution.hess);
+
+        using MatrixType = typename HPolytope::MT;
+        using Input = crhmc_input
+            <
+                MatrixType,
+                Point,
+                NegativeLogprobFunctor,
+                NegativeGradientFunctor,
+                HessianFunctor
+            >;
+        Input input = convert2crhmc_input<Input>(HP, F, G, H);
+
+        using CrhmcProblem = crhmc_problem<Point, Input>;
+        CrhmcProblem problem = CrhmcProblem(input);
+        if(problem.terminate){return;}
+
+        using Solver = ImplicitMidpointODESolver
+            <
+                Point, NT, CrhmcProblem,
+                NegativeGradientFunctor, simdLen
+            >;
+
+        using Walk = typename WalkType::template Walk
+            <
+                    Point,
+                    CrhmcProblem,
+                    RandomNumberGenerator,
+                    NegativeGradientFunctor,
+                    NegativeLogprobFunctor,
+                    Solver
+            >;
+        using WalkParams = typename WalkType::template parameters
+            <
+                    NT,
+                    NegativeGradientFunctor
+            >;
+        Point p = Point(problem.center);
+        problem.options.simdLen = simdLen;
+        WalkParams params(distribution.L, p.dimension(), problem.options);
+
+        // TODO: pass a user defined eta to bypass the one created by the walk using L
+        //if (distribution.eta > 0) {
+        //    params.eta = distribution.eta;
+        // }
+
+        Walk walk(problem, p, input.df, input.f, params);
+
+        for (unsigned int i = 0; i < nburns; ++i)
+        {
+            walk.apply(rng, walk_len);
+        }
+
+        for (unsigned int i = 0; i < std::ceil((float)rnum/simdLen); ++i)
+        {
+            walk.apply(rng, walk_len);
+            if (walk.P.terminate) {return;}
+            auto x = walk.getPoints();
+
+            if ((i + 1) * simdLen > rnum)
+            {
+                for (int j = 0; j < rnum-simdLen*i; j++)
+                {
+                    samples.col(i+j) = x.col(j);
+                }
+                break;
+            }
+            for (int j = 0; j < x.cols(); j++)
+            {
+                samples.col(i+j) = x.col(j);
+            }
+        }
+    }
+
+    else
+    {
+        static_assert(detail::AlwaysFalse<WalkType, Distribution>::value,
+            "Not implemented for this combination of walk and distribution.");
+    }
+}
+
+
+#endif //SAMPLE_POINTS_HPP
\ No newline at end of file