diff --git a/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h
index b4a69ce98..809994c04 100644
--- a/include/convex_bodies/hpolytope.h
+++ b/include/convex_bodies/hpolytope.h
@@ -7,6 +7,7 @@
 //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.
+//Contributed and/or modified by Iva Janković, as part of Google Summer of Code 2025 program.
 
 // Licensed under GNU LGPL.3, see LICENCE file
 
@@ -203,6 +204,13 @@ class HPolytope {
         return b;
     }
 
+    VT get_row(int facet_index) const 
+    {
+        if (facet_index < 0 || facet_index >= A.rows())
+            throw std::out_of_range("Facet index out of range!");
+        return A.row(facet_index);
+    }
+    
     bool is_normalized ()
     {
         return normalized;
@@ -535,6 +543,7 @@ class HPolytope {
         NT inner_prev = params.inner_vi_ak;
         VT sum_nom;
         int m = num_of_hyperplanes(), facet;
+        int skip = params.facet_prev;
 
         Ar.noalias() += lambda_prev*Av;
         if(params.hit_ball) {
@@ -547,20 +556,16 @@ class HPolytope {
         NT* sum_nom_data = sum_nom.data();
         NT* Av_data = Av.data();
 
-        for (int i = 0; i < m; i++) {
-            if (*Av_data == NT(0)) {
-                //std::cout<<"div0"<<std::endl;
-                ;
-            } else {
-                lamda = *sum_nom_data / *Av_data;
-                if (lamda < min_plus && lamda > 0) {
-                    min_plus = lamda;
-                    facet = i;
-                    params.inner_vi_ak = *Av_data;
-                }
+        for (int i = 0; i < m; ++i, ++Av_data, ++sum_nom_data){
+            if (i == skip) continue;    
+            if (*Av_data == NT(0)) continue;
+
+            lamda = *sum_nom_data / *Av_data;
+            if (lamda < min_plus && lamda > 0) {
+                min_plus = lamda;
+                facet = i;
+                params.inner_vi_ak = *Av_data;
             }
-            Av_data++;
-            sum_nom_data++;
         }
         params.facet_prev = facet;
         return std::pair<NT, int>(min_plus, facet);
@@ -641,6 +646,7 @@ class HPolytope {
         NT lamda = 0;
         VT sum_nom;
         int m = num_of_hyperplanes(), facet;
+        int skip = params.facet_prev;
 
         Ar.noalias() += lambda_prev*Av;
         sum_nom.noalias() = b - Ar;
@@ -649,24 +655,20 @@ class HPolytope {
         NT* sum_nom_data = sum_nom.data();
         NT* Av_data = Av.data();
 
-        for (int i = 0; i < m; i++) {
-            if (*Av_data == NT(0)) {
-                //std::cout<<"div0"<<std::endl;
-                ;
-            } else {
-                lamda = *sum_nom_data / *Av_data;
-                if (lamda < min_plus && lamda > 0) {
-                    min_plus = lamda;
-                    facet = i;
-                    params.inner_vi_ak = *Av_data;
-                }
+        for (int i = 0; i < m; ++i, ++Av_data, ++sum_nom_data) {
+            if (i == skip) continue;
+            if (*Av_data == NT(0)) continue;
+            lamda = *sum_nom_data / *Av_data;
+            if (lamda < min_plus && lamda > 0) {
+                min_plus = lamda;
+                facet = i;
+                params.inner_vi_ak = *Av_data;
             }
-            Av_data++;
-            sum_nom_data++;
         }
         params.facet_prev = facet;
         return std::pair<NT, int>(min_plus, facet);
     }
+    
 
     //-----------------------------------------------------------------------------------//
 
@@ -1086,4 +1088,4 @@ class HPolytope {
     }
 };
 
-#endif
\ No newline at end of file
+#endif
diff --git a/include/diagnostics/scaling_ratio.hpp b/include/diagnostics/scaling_ratio.hpp
new file mode 100644
index 000000000..b97953677
--- /dev/null
+++ b/include/diagnostics/scaling_ratio.hpp
@@ -0,0 +1,126 @@
+// VolEsti (volume computation and sampling library)
+
+// Copyright (c) 2012-2025 Vissarion Fisikopoulos
+// Copyright (c) 2018-2025 Apostolos Chalkis
+// Copyright (c) 2025-2025 Iva Janković
+
+// Contributed and/or modified by Iva Janković, as part of Google Summer of Code 2025 program.
+
+// Licensed under GNU LGPL.3, see LICENCE file
+
+#ifndef DIAGNOSTICS_SCALING_RATIO_HPP
+#define DIAGNOSTICS_SCALING_RATIO_HPP
+
+#include <tuple>
+#include <vector>
+#include <cmath>
+#include <cstdlib>
+
+template<typename Polytope>
+std::tuple<typename Polytope::VT, typename Polytope::MT,typename Polytope::VT, typename Polytope::VT>
+scaling_ratio_boundary_test(const Polytope&  P,
+                            const typename Polytope::MT& samples,
+                            const typename Polytope::NT tol = 1e-10,
+                            const typename Polytope::NT min_ratio = 0.01) 
+{
+    using VT = typename Polytope::VT;
+    using MT = typename Polytope::MT;
+    using NT = typename Polytope::NT;
+    
+    const int dim = P.dimension();
+    const int m = P.num_of_hyperplanes();
+    const unsigned n_samp = static_cast<unsigned>(samples.cols());
+
+    VT scale(10); //we scale 10 times
+    MT coverage(m, 10);
+    std::vector<int> facet_id(n_samp, -1);
+    const auto A_full = P.get_mat();
+    const auto b_full = P.get_vec();
+
+    for (int i = 0; i < n_samp; ++i) {
+        auto Aq = A_full * samples.col(i);
+        
+        for (size_t k = 0; k < m; ++k) {
+            if (std::abs(Aq[k] - b_full[k]) < tol) 
+            {
+                facet_id[i] = static_cast<int>(k);
+                break;
+            }
+        }
+    }
+
+    for (int f = 0; f < m; ++f) 
+    {
+        // Samples S on facet f
+        std::vector<int> S;
+        S.reserve(n_samp / m);
+        for (unsigned i = 0; i < n_samp; ++i) {
+            if (facet_id[i] == f)
+                S.push_back(static_cast<int>(i));
+        }
+
+        const double ratio = static_cast<double>(S.size()) / n_samp;
+        if (ratio < min_ratio)  continue;
+
+        // Finding the center
+        VT p = VT::Zero(dim);
+        for (int idx : S) p += samples.col(idx);
+        p /= static_cast<double>(S.size());
+
+        // Looping over scale factor
+        for (int k = 0; k < 10; ++k) 
+        {
+            NT step = 0.1 * (k+1);
+            NT x = std::pow(step, 1.0 / dim);
+            // Local copy of polytope for each scaling
+            Polytope P_loc = P;
+            
+            //Shifting and scaling
+            P_loc.shift(p);
+            MT T = (1.0 / x) * MT::Identity(dim, dim);
+            P_loc.linear_transformIt(T);
+
+            //Parameters of new polytope
+            const auto& A_sh = P_loc.get_mat();
+            const auto& b_sh = P_loc.get_vec();
+
+            // Points still in facet
+            unsigned survivors = 0;
+            for (int idx : S) {
+
+                const VT q_shift = samples.col(idx) - p;
+                bool inside = true;
+
+                for (int j = 0; j < A_sh.rows(); ++j) {
+                    if (j == f) continue;
+                    if (A_sh.row(j).dot(q_shift) - b_sh[j] > tol) {
+                        inside = false; break;
+                    }
+                }
+                if (inside) ++survivors;
+            }
+
+            coverage(f, k) = double(survivors) / double(S.size());
+            scale[k]=std::pow(x, dim);
+        }
+    }
+
+    int K = scale.size();           // number of scaling factors
+    VT max_dev(m), avg_dev(m);      // Vectors of average and maximum deviations
+
+    for (int f = 0; f < m; ++f) {
+        double sumd = 0.0;
+        double maxd = 0.0;
+        for (int k = 0; k < K; ++k) {
+            double d = std::abs(coverage(f, k) - scale[k])* 100.0; // in percentage
+            sumd += d;
+            if (d > maxd) maxd = d;
+        }
+        avg_dev[f] = sumd/K;
+        max_dev[f] = maxd;      
+    }
+
+    return {scale, coverage, max_dev, avg_dev};
+}
+
+#endif // DIAGNOSTICS_SCALING_RATIO_HPP
diff --git a/include/preprocess/feasible_point.hpp b/include/preprocess/feasible_point.hpp
index 6e7f5db0a..c151586f2 100644
--- a/include/preprocess/feasible_point.hpp
+++ b/include/preprocess/feasible_point.hpp
@@ -30,5 +30,41 @@ VT compute_feasible_point(MT const& A, VT const& b)
     return x;
 }
 
+template <typename Point, typename Polytope, typename RNG>
+std::pair<typename Polytope::VT,int> compute_feasible_boundary_point(Polytope const& P,
+                                                                     RNG& rng,
+                                                                     typename Point::FT eps)
+{
+    using VT = typename Polytope::VT;
+    using NT = typename Point::FT;
+    const std::size_t m = P.num_of_hyperplanes();
+
+    // Find the interior point
+    VT r = compute_feasible_point(P.get_mat(), P.get_vec());
+
+    // Random ray 
+    const int dim = P.dimension();
+    Point v_pt = GetDirection<Point>::apply(dim, rng);
+    VT v = v_pt.getCoefficients();
+
+    // First‐hit oracle
+    VT Ar(m), Av(m);
+    struct Params { NT inner_vi_ak; int facet_prev; };
+    Params params;
+
+    auto [lambda_min, facet] = P.line_first_positive_intersect(r, v, Ar, Av, params);
+
+    // Checks
+    if (!std::isfinite(lambda_min) || lambda_min <= eps || facet < 0)
+        throw std::runtime_error("Failed to hit boundary!!!");
+
+    // Compute boundry point + final check
+    VT x = r + lambda_min * v;
+    if ((P.get_mat() * x - P.get_vec()).maxCoeff() > eps)
+        throw std::runtime_error("Boundary point violates constraints!!!");
+
+    return std::pair<VT,int>(x, facet);    
+}
+
 
 #endif
diff --git a/include/random_walks/shake_and_bake_walk.hpp b/include/random_walks/shake_and_bake_walk.hpp
new file mode 100644
index 000000000..58267240e
--- /dev/null
+++ b/include/random_walks/shake_and_bake_walk.hpp
@@ -0,0 +1,167 @@
+// VolEsti (volume computation and sampling library)
+
+// Copyright (c) 2012-2025 Vissarion Fisikopoulos
+// Copyright (c) 2018-2025 Apostolos Chalkis
+// Copyright (c) 2025-2025 Iva Janković
+
+// Contributed and/or modified by Iva Janković, as part of Google Summer of Code 2025 program.
+
+// Licensed under GNU LGPL.3, see LICENCE file
+
+#ifndef RANDOM_WALKS_SHAKE_AND_BAKE_WALK_HPP
+#define RANDOM_WALKS_SHAKE_AND_BAKE_WALK_HPP
+
+#include <Eigen/Eigen>
+#include <Eigen/Dense>
+#include <cmath>
+#include <algorithm>
+#include <stdexcept> 
+
+#include "sampling/sphere.hpp"
+#include "convex_bodies/hpolytope.h"
+#include "convex_bodies/correlation_matrices/corre_matrix.hpp"
+
+struct ShakeAndBakeWalk
+{
+
+    template <typename Polytope, typename RandomNumberGenerator>
+    struct Walk
+    {
+        using Point = typename Polytope::PointType;
+        using VT = typename Polytope::VT;
+        using NT = typename Point::FT;
+        using MT = typename Polytope::MT;
+
+        struct update_parameters 
+        {
+            int facet_prev = -1;
+            NT  inner_vi_ak = NT(0);   
+        };
+
+        update_parameters _params;
+
+        static constexpr NT kDefaultEpsilon = NT(1e-10);
+
+        template <typename GenericPolytope>
+        Walk(GenericPolytope& P,
+            const Point& boundary_pt, 
+            int facet_idx,    
+            RandomNumberGenerator& rng,
+            NT eps = kDefaultEpsilon)
+            : epsilon_{eps}
+        {
+            P.normalize();
+            initialize(P,boundary_pt, facet_idx, rng);
+        }
+
+        NT   get_epsilon() const noexcept { return epsilon_; }
+
+        void apply(Polytope& P, unsigned int walk_len, RandomNumberGenerator& rng)
+        {
+            const NT eps = epsilon_; 
+
+            for (unsigned step = 0; step < walk_len; ++step)
+            {
+                Point v = get_direction(P,rng);
+
+                int facet_new;
+                std::tie(_lambda_hit, facet_new) = P.line_positive_intersect(_p, v, _Ar, _Av, _lambda_hit, _params);
+
+                if (!std::isfinite(_lambda_hit) || _lambda_hit <= eps  || facet_new < 0) 
+                {
+                    _lambda_hit = NT(0);
+                    continue;
+                }
+
+                _p += _lambda_hit * v;
+                _A_row_k = P.get_row(facet_new);
+                _params.facet_prev  = facet_new;
+            }
+        }
+
+
+        const Point& getCurrentPoint() const noexcept { return _p; }
+
+    private:
+
+        Point get_direction(Polytope& P, RandomNumberGenerator& rng)
+        {
+            int _dim = P.dimension();
+            VT z = GetDirection<Point>::apply(_dim, rng).getCoefficients();
+            MT I_cc = - _A_row_k * _A_row_k.transpose();
+            I_cc.diagonal() += VT::Ones(_dim);
+            NT U = rng.sample_urdist();               
+            NT r = std::pow(U, NT(1)/NT(_dim-1)); 
+            NT cz = _A_row_k.dot(z);
+            VT z_tilde  = I_cc*z;
+            z_tilde *= r;
+            z_tilde /= std::sqrt(NT(1) - cz*cz);
+            
+            VT v = z_tilde - std::sqrt(NT(1) - r*r) * _A_row_k;
+            return Point(v);
+        }
+
+        void initialize(Polytope& P,
+                        const Point& boundary_pt,
+                        int facet_idx,
+                        RandomNumberGenerator& rng)
+        {
+            int _dim = P.dimension();
+            int _m = P.num_of_hyperplanes();
+            VT b=P.get_vec();
+            int active_facet;
+
+            NT kFacetEps = epsilon_;
+
+            // Checking if boundary point belongs to facet_idx
+            _p = boundary_pt;
+            VT ai = P.get_row(facet_idx);
+            NT dist = std::abs(ai.dot(_p.getCoefficients()) - b.coeff(facet_idx));
+            if (dist > kFacetEps)
+            {
+                active_facet = -1;
+                for (int i = 0; i < _m; ++i) 
+                {
+                    VT ai = P.get_row(i);
+                    NT dist = std::abs(ai.dot(_p.getCoefficients()) - b.coeff(i));
+                    if (dist < kFacetEps) 
+                    {
+                        active_facet = i;
+                        break;
+                    }
+                }
+                if (active_facet < 0)
+                {
+                    throw std::runtime_error("Boundary point not on any facet!");
+                }
+            }
+            active_facet = facet_idx;
+
+            //Normal of active facet
+            _A_row_k = P.get_row(active_facet);
+
+            //Initializing Ar and Av 
+            _Ar.setZero(_m);
+            _Av.setZero(_m);
+            _lambda_hit = NT(0);
+            
+            //Calculating first Ar and initializing lambda 
+            _Ar.noalias() = P.get_mat() * _p.getCoefficients();
+            _lambda_hit = NT(0);
+
+            _A_row_k = P.get_row(active_facet);
+
+            _params.facet_prev = active_facet;
+
+        }
+
+        NT epsilon_{kDefaultEpsilon};
+        Point _p;
+        VT _Ar;            
+        VT _Av;            
+        NT _lambda_hit;
+        VT _A_row_k;
+    };
+};
+
+#endif // RANDOM_WALKS_SHAKE_AND_BAKE_WALK_HPP
diff --git a/include/sampling/random_point_generators.hpp b/include/sampling/random_point_generators.hpp
index 998d8dafc..f15bae891 100644
--- a/include/sampling/random_point_generators.hpp
+++ b/include/sampling/random_point_generators.hpp
@@ -222,6 +222,33 @@ struct BoundaryRandomPointGenerator
     }
 };
 
+template <typename Walk>
+struct ShakeAndBakeRandomPointGenerator
+{
+    template <
+        typename Polytope,
+        typename Point, 
+        typename PointList,
+        typename WalkPolicy,
+        typename RandomNumberGenerator
+    >
+    static void apply(Polytope& P,
+                      Point &p,
+                      unsigned int rnum,
+                      unsigned int walk_len,
+                      PointList& randPoints,
+                      WalkPolicy& policy,
+                      RandomNumberGenerator& rng,
+                      int facet_idx = -1) 
+    {
+        Walk walk(P, p, facet_idx, rng);
+
+        for (unsigned int i = 0; i < rnum; ++i) {
+            walk.apply(P, walk_len, rng);
+            policy.apply(randPoints, walk.getCurrentPoint());
+        }
+    }
+};
 
 template
 <
diff --git a/include/sampling/sampling.hpp b/include/sampling/sampling.hpp
index 696186098..8c789bb33 100644
--- a/include/sampling/sampling.hpp
+++ b/include/sampling/sampling.hpp
@@ -220,6 +220,39 @@ void gaussian_sampling(PointList &randPoints,
                                 push_back_policy, rng, WalkType.param);
 }
 
+template <
+    typename WalkTypePolicy,
+    typename PointList,
+    typename Polytope,
+    typename RandomNumberGenerator,
+    typename Point
+>
+void shakeandbake_sampling(PointList &randPoints,
+                           Polytope &P,
+                           RandomNumberGenerator &rng,
+                           const unsigned int walk_len,
+                           const unsigned int  rnum,
+                           const Point &starting_point,
+                           unsigned int const  &nburns = 0,
+                           int facet_idx = -1)  
+{
+    typedef typename WalkTypePolicy::template Walk
+            <
+                    Polytope,
+                    RandomNumberGenerator
+            > walk;
+
+    PushBackWalkPolicy push_policy;
+
+    if (nburns > 0) 
+    {
+        ShakeAndBakeRandomPointGenerator<walk>::apply(P,starting_point,nburns,walk_len,randPoints,push_policy,rng,facet_idx);     
+        randPoints.clear();    
+    }
+
+    ShakeAndBakeRandomPointGenerator<walk>::apply(P,starting_point,rnum,walk_len,randPoints,push_policy,rng,facet_idx);
+}
+
 template <
         typename PointList,
         typename Polytope,
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 85c5fa56d..58b12a693 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -212,6 +212,10 @@ add_test(NAME test_ghmc COMMAND sampling_test -tc=ghmc)
 add_test(NAME test_gabw COMMAND sampling_test -tc=gabw)
 add_test(NAME test_sparse COMMAND sampling_test -tc=sparse)
 
+add_executable (shake_and_bake_test shake_and_bake_test.cpp $<TARGET_OBJECTS:test_main>)
+add_test(NAME test_shake_and_bake COMMAND shake_and_bake_test -tc=shake_and_bake)
+add_test(NAME test_scaling_ratio_boundary COMMAND shake_and_bake_test -tc=scaling_ratio_boundary)
+
 add_executable (mmcs_test mmcs_test.cpp $<TARGET_OBJECTS:test_main>)
 add_test(NAME test_mmcs COMMAND mmcs_test -tc=mmcs)
 
@@ -404,6 +408,7 @@ TARGET_LINK_LIBRARIES(volume_cb_vpoly_intersection_vpoly lp_solve ${MKL_LINK} co
 TARGET_LINK_LIBRARIES(rounding_test lp_solve ${MKL_LINK} coverage_config)
 TARGET_LINK_LIBRARIES(mcmc_diagnostics_test lp_solve ${MKL_LINK} coverage_config)
 TARGET_LINK_LIBRARIES(sampling_test lp_solve ${MKL_LINK} coverage_config)
+TARGET_LINK_LIBRARIES(shake_and_bake_test lp_solve ${MKL_LINK} coverage_config)
 TARGET_LINK_LIBRARIES(mmcs_test lp_solve ${MKL_LINK} coverage_config)
 TARGET_LINK_LIBRARIES(benchmarks_sob lp_solve ${MKL_LINK} coverage_config)
 TARGET_LINK_LIBRARIES(benchmarks_cg lp_solve ${MKL_LINK} coverage_config)
diff --git a/test/shake_and_bake_test.cpp b/test/shake_and_bake_test.cpp
new file mode 100644
index 000000000..bfc38f0d0
--- /dev/null
+++ b/test/shake_and_bake_test.cpp
@@ -0,0 +1,133 @@
+// VolEsti (volume computation and sampling library)
+
+// Copyright (c) 2012-2025 Vissarion Fisikopoulos
+// Copyright (c) 2018-2025 Apostolos Chalkis
+// Copyright (c) 2025-2025 Iva Janković
+
+// Contributed and/or modified by Iva Janković, as part of Google Summer of Code 2025 program.
+
+// Licensed under GNU LGPL.3, see LICENCE file
+
+#include "doctest.h"
+#include <fstream>
+#include <iostream>
+
+#include <boost/random.hpp>
+#include <boost/random/uniform_int.hpp>
+#include <boost/random/normal_distribution.hpp>
+#include <boost/random/uniform_real_distribution.hpp>
+
+#include "misc/misc.h"
+#include "random_walks/shake_and_bake_walk.hpp"
+
+#include "volume/volume_sequence_of_balls.hpp"
+#include "generators/known_polytope_generators.h"
+#include "sampling/sampling.hpp"
+
+#include "diagnostics/univariate_psrf.hpp"
+#include "diagnostics/scaling_ratio.hpp"
+
+#include "preprocess/feasible_point.hpp"
+
+template
+<
+    typename MT,
+    typename WalkType,
+    typename Polytope
+>
+
+MT get_samples_shake_and_bake(Polytope &P)
+{
+    typedef typename Polytope::PointType Point;
+    typedef typename Polytope::NT NT;
+
+    typedef BoostRandomNumberGenerator<boost::mt19937, NT, 3> RNGType;
+
+    unsigned int walkL = 10, numpoints = 10000, nburns = 0, d = P.dimension();
+    RNGType rng(d);
+
+    auto [boundary_pt, facet_idx] = compute_boundary_point<Point>(P, rng, 1e-7);
+
+    std::list<Point> randPoints;
+
+    shakeandbake_sampling<WalkType>(randPoints,P,rng,walkL,numpoints,boundary_pt, nburns,facet_idx);        
+
+    MT samples(d, numpoints);
+
+    unsigned int jj = 0;
+    for (const Point& q : randPoints)
+    {
+        samples.col(jj++) = q.getCoefficients();
+    }
+        
+    return samples;
+}
+
+
+template <typename NT, typename WalkType = ShakeAndBakeWalk>
+void call_test_shake_and_bake(){
+    typedef Cartesian<NT>    Kernel;
+    typedef typename Kernel::Point    Point;
+    typedef HPolytope<Point> Hpolytope;
+    typedef Eigen::Matrix<NT,Eigen::Dynamic,Eigen::Dynamic> MT;
+    typedef Eigen::Matrix<NT,Eigen::Dynamic,1> VT;
+    Hpolytope P;
+    unsigned int d = 10;
+
+    std::cout << "--- Testing Running Shake and Bake for H-cube 10" << std::endl;
+    P = generate_cube<Hpolytope>(d, false);
+    P.ComputeInnerBall();
+
+    MT samples = get_samples_shake_and_bake<MT, WalkType>(P);
+
+    VT score = univariate_psrf<NT, VT>(samples);
+    std::cout << "psrf = " << score.maxCoeff() << std::endl;
+
+    CHECK(score.maxCoeff() < 1.1);
+}
+
+template <typename NT, typename WalkType = ShakeAndBakeWalk>
+void call_test_scaling_ratio_boundary(){
+    typedef Cartesian<NT>    Kernel;
+    typedef typename Kernel::Point    Point;
+    typedef HPolytope<Point> Hpolytope;
+    typedef Eigen::Matrix<NT,Eigen::Dynamic,Eigen::Dynamic> MT;
+    typedef Eigen::Matrix<NT,Eigen::Dynamic,1> VT;
+    Hpolytope P;
+    unsigned int d = 10;
+
+    std::cout << "--- Testing Scaling Ratio Boundary Test for H-cube 10" << std::endl;
+    P = generate_cube<Hpolytope>(d, false);
+    P.ComputeInnerBall();
+
+    // Get boundary samples using shake and bake
+    MT samples = get_samples_shake_and_bake<MT, WalkType>(P);
+
+    // Call scaling_ratio_boundary_test
+    auto [scale, coverage, max_dev, avg_dev] = scaling_ratio_boundary_test(P, samples);
+    
+    std::cout << "Scale factors: ";
+    for (int i = 0; i < scale.size(); ++i) {
+        std::cout << scale[i] << " ";
+    }
+    std::cout << std::endl;
+    
+    // Check --for each facet-- that each value in coverage is at most 0.1 different from the corresponding value in scale
+    const int m = P.num_of_hyperplanes();
+    const int K = scale.size();
+    
+    for (int f = 0; f < m; ++f) {
+        for (int k = 0; k < K; ++k) {
+            NT diff = std::abs(coverage(f, k) - scale[k]);
+            CHECK(diff <= 0.01);
+        }
+    }
+}
+
+TEST_CASE("shake_and_bake") {
+    call_test_shake_and_bake<double>();
+}
+
+TEST_CASE("scaling_ratio_boundary") {
+    call_test_scaling_ratio_boundary<double>();
+}