Shake and bake sampling (GeomScale#361)

* Final reorganization Shake and Bake

Signed-off-by: Iva Janković <ivaxjankovic@gmail.com>

* Update shakeandbake.cpp with improved diagnostics

* Updates for high dimensional polytopes

* Uniformity tests updates

* For debugging: Sampling testing

* Corrected testing

* Updated shake and bake walk with some corrections

* Apply func fixes

* Apply func fixes in hpp

* Compute boundary point function

* Revised shake and bake

* Just running SB

* Everything for revision

* Quick fix of SB

* 2 helper functions

* Fixed SB

* Line_positive_intersect_skip func

* For revision

* Latest version SB

* SB without v*=1

* Latest hpolytope.h

* Get_direction romeijn1991

* SB fixed + scaling_ratio.hpp

* SB + scaling test fixed

* Fixed scaling_ratio.hpp

* Latest SB (hpp and cpp)

* Ready to merge sb.hpp?

* For revision 2

* For Revision 2.1

* For revision 3

* Scaling ratio test fixed

* Cleaning 1

* Cleaned 2

* Ready to go!

* Shake and bake: unit tests

* Unit test: polished!

* Update README.md

---------

Signed-off-by: Iva Janković <ivaxjankovic@gmail.com>

Author: ivajankovic

examples/sb_sampling/CMakeLists.txt

@@ -0,0 +1,103 @@
+# 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
+
+project( VolEsti )
+
+
+CMAKE_MINIMUM_REQUIRED(VERSION 3.11)
+
+set(CMAKE_ALLOW_LOOSE_LOOP_CONSTRUCTS true)
+
+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)
+
+
+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()
+
+
+# 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")
+
+  include_directories (BEFORE ../../external)
+  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
+  add_definitions(${CMAKE_CXX_FLAGS} "-O3")  # optimization of the compiler
+  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 (shakeandbake shakeandbake.cpp)
+  TARGET_LINK_LIBRARIES(shakeandbake ${LP_SOLVE})
+
+endif()

examples/sb_sampling/README.md

@@ -0,0 +1,31 @@
+## Explanation
+In 'shake_and_bake_walk.hpp' the Running variant of Shake And Bake class of boundary sampling algorithms. It follows the steps as described in [1]. The walk can be run by using 'shakeandbake.cpp' Additionaly, to test the uniformity, the scaling ratio test ('scaling_ratio.hpp') along with random facet 2D projection for generic polytopes (cube, simplex, birkhoff) with the plots are added in 'sbtest.py' (with generated txt files). 
+
+[1] C. G. E. Boender, R. J. Caron, J. F. McDonald, A. H. G. Rinnooy Kan,  
+    H. E. Romeijn, R. L. Smith, J. Telgen i A. C. F. Vorst,  
+    *Shake-And-Bake Algorithms for Generating Uniform Points on the Boundary of Bounded Polyhedra*, 1991.  
+    Available at: https://doi.org/10.1016/0166-218X(91)90006-7
+
+## Original and Limping Variant
+The implemented Shake And Bake variant is Running SB because it performs significantly better than its counterparts Original and Limping. But, if one wants to test that out, here is the additional piece of code to be added in `apply` function alongside some simple enum switch / branch logic. Also, it is very important to use  `_v = GetDirection<Point>::apply(n, rng); ` (coressponds to Step 1. of Original and Limping SB from the paper ) instead of `Point v = get_direction(P,rng);` (corresponds to Step 1. for Running SB).  
+
+```bash
+NT beta;
+if (mode_ == Mode::Original) {
+    NT den = dot_r - dot_k;
+    if (std::abs(den) < eps) continue;
+    beta = std::clamp(dot_r / den, NT(0), NT(1));
+} 
+else {
+    beta = -dot_k;
+}
+
+if (beta > NT(0) && beta <= NT(1) &&
+    rng.sample_urdist() < beta)
+{
+    p_         = y;
+    facet_idx_ = facet_new;
+    A_row_k_   = A_row_r;
+    Ar_.noalias() -= lambda_hit * Av_;  
+}
+```  

examples/sb_sampling/sbtest.py

@@ -0,0 +1,97 @@
+# 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
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+shape = input("Name of the polytope ('cube', 'simplex', 'birkhoff', or anything else): ").strip()
+dim   = int(input("Number of dimensions: ").strip())
+
+# Loading the data
+filename = f"build/sb_{shape}_{dim}_run.txt"
+data     = np.loadtxt(filename)
+n, d     = data.shape
+eps = 1e-7
+
+# Finding the facet
+if shape == "cube":
+    facet_col = np.random.randint(0, d)
+    facet_val = np.random.choice([+1.0, -1.0])
+    desc      = f"x_{facet_col}≈{facet_val}"
+    mask      = np.abs(data[:, facet_col] - facet_val) < eps
+
+elif shape in ("birkhoff", "simplex"):
+    facet_col = np.random.randint(0, d)
+    facet_val = 0.0
+    desc      = f"x_{facet_col}≈0"
+    mask      = np.abs(data[:, facet_col] - facet_val) < eps
+
+else:
+    facet_col = None
+    desc      = "full cloud"
+    mask      = np.ones(n, dtype=bool)
+
+facet_pts = data[mask]
+
+# "Center" of the facet
+centroid = facet_pts.mean(axis=0)
+
+# Projection on two random axes
+axes = [i for i in range(d) if i != facet_col] or list(range(d))
+i, j = np.random.choice(axes, size=2, replace=False)
+
+# Scatter of points in the facet
+
+plt.figure(figsize=(6,6))
+plt.scatter(facet_pts[:, i], facet_pts[:, j], s=12, alpha=0.7, label="Samples")
+plt.scatter(centroid[i], centroid[j], c='red', marker='*', s=150, label='Center')
+plt.xlabel(f"x_{i}")
+plt.ylabel(f"x_{j}")
+plt.title(f"{shape.capitalize()} {dim}D — facet {desc} with {len(facet_pts)} points ")
+plt.legend()
+plt.grid(True)
+plt.tight_layout()
+plt.show()
+
+# Scaling ratio testing
+coverage = {}
+filename = f"build/sb_{shape}_{dim}_coverage.txt"
+with open(filename) as f:
+    for line in f:
+        line = line.strip()
+        if not line:
+            continue
+        parts = line.split()
+        if parts[0] != 'Facet':
+            continue
+        facet = int(parts[1])
+        tail = ' '.join(parts[4:])
+        pairs = [p.strip().rstrip(',') for p in tail.split(',') if p.strip()]
+        xs, covs = [], []
+        for p in pairs:
+            x_str, cov_str = p.split(':', 1)
+            x = float(x_str)
+            cov = float(cov_str)
+            xs.append(x)
+            covs.append(cov)
+        coverage[facet] = (np.array(xs), np.array(covs))
+
+plt.figure(figsize=(8,6))
+plt.grid(True) 
+for facet, (xs, covs) in coverage.items():
+    plt.plot(xs, covs, label=f'Facet {facet}')
+
+plt.xlabel(r'$x^{\,d}$')
+plt.ylabel('Coverage ratio')
+plt.grid(True) 
+plt.title('Scaling coverage per facet (vs $x^{d}$)')
+plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
+plt.tight_layout()
+plt.show()

examples/sb_sampling/shakeandbake.cpp

@@ -0,0 +1,159 @@
+// 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 <iostream>
+#include <fstream>
+#include <string>
+#include <vector>
+#include <Eigen/Eigen>
+#include <boost/random.hpp>
+
+#include "cartesian_geom/cartesian_kernel.h"
+#include "preprocess/feasible_point.hpp"   
+#include "convex_bodies/hpolytope.h"
+#include "random_walks/random_walks.hpp"
+#include "random_walks/shake_and_bake_walk.hpp"
+#include "generators/known_polytope_generators.h"
+#include "diagnostics/scaling_ratio.hpp" 
+
+
+using NT = double;
+using Kernel = Cartesian<NT>;
+using Point = Kernel::Point;
+using RNG  = BoostRandomNumberGenerator<boost::random::mt19937, NT>;
+using HPoly = HPolytope<Point>;
+using Walker1 = ShakeAndBakeWalk::Walk<HPoly, RNG>;
+
+
+int main(int argc, char* argv[])
+{
+    if (argc < 3) {
+        std::cerr << "Usage: " << argv[0]
+                  << " <cube|simplex|birkhoff> <dimension> [epsilon]\n";
+        return 1;
+    }
+
+    std::string shape  = argv[1];
+    unsigned  cli_n = std::stoi(argv[2]);
+    NT eps_cli = (argc > 3)
+                 ? static_cast<NT>(std::stod(argv[3]))
+                 : Walker1::kDefaultEpsilon;      // default value if not manually
+
+    //Generating polytope 
+    HPoly P;
+    if (shape == "cube") P = generate_cube<HPoly>(cli_n, false);
+    else if (shape == "simplex") P = generate_simplex<HPoly>(cli_n, false);
+    else if (shape == "birkhoff") P = generate_birkhoff<HPoly>(cli_n);
+    else {
+        std::cerr << "Unknown polytope type: " << shape << '\n';
+        return 1;
+    }
+
+    // Walk parameters adjustments
+    const unsigned true_dim = P.dimension();
+    unsigned walk_len, n_samples, burn_in_iters;
+
+    int mode = (shape == "cube" || shape == "simplex") ? 0
+         : (shape == "birkhoff")             ? 1
+         : 2;
+
+    switch (mode) {
+        case 0:  // cube or simplex
+            walk_len  = 20 * true_dim;
+            n_samples = 500 * true_dim;
+            burn_in_iters = 5  * true_dim;
+            break;
+
+        case 1:  // birkhoff
+            walk_len = 100 * true_dim;
+            n_samples  = 2000 * true_dim;
+            burn_in_iters = 10  * true_dim;
+            break;
+
+        default: 
+            walk_len = 20 * true_dim;
+            n_samples  = 100 * true_dim;
+            burn_in_iters = 20 * true_dim;
+            break;
+    }
+
+    std::cout << "Parameters: walk_len="   << walk_len
+              << ", n_samples="            << n_samples
+              << ", burn_in_iters="        << burn_in_iters
+              << " (dim="                  << true_dim
+              << ") eps="                  << eps_cli << '\n';
+
+    // Initializing the walk
+    RNG rng(true_dim);
+
+    auto [boundary_pt, facet_idx] = compute_boundary_point<Point>(P, rng, eps_cli);
+    Walker1 walk1(P, boundary_pt, facet_idx, rng,eps_cli);
+    const NT tol = walk1.get_epsilon();                     
+
+    const std::string base = "sb_" + shape + "_" + std::to_string(cli_n);
+    std::ofstream out(base + "_run.txt");
+
+
+    Eigen::Matrix<NT, Eigen::Dynamic, Eigen::Dynamic> samples1(true_dim, n_samples);
+    std::vector<int> facet_id1(n_samples, -1);
+
+    //Burn in 
+    for (int i = 0; i < burn_in_iters; ++i)
+        walk1.apply(P, walk_len, rng);
+
+    // Sampling
+    for (int i = 0; i < n_samples; ++i) {
+        walk1.apply(P, walk_len, rng);
+        const Point& q = walk1.getCurrentPoint();
+        samples1.col(i) = q.getCoefficients();
+
+        // File inscription
+        for (unsigned d = 0; d < true_dim; ++d)
+            out << q[d] << (d + 1 < true_dim ? ' ' : '\n');
+    }
+    out.close();
+    
+    std::cout << "Generated " << n_samples << " samples in "
+            << walk_len   << " steps each.\n";
+
+
+    //Scaling ratio test 
+    auto [scales, coverage, max_dev, avg_dev] = scaling_ratio_boundary_test(P, samples1,tol);
+    
+    std::cout << "Scaling factors:\n";
+    for (double s : scales) {
+        std::cout << s << " ";
+    }
+    std::cout << "\n\nCoverage matrix (each row = one facet):\n";
+    for (int f = 0; f < coverage.rows(); ++f) {
+        std::cout << "Facet " << f << ": ";
+        for (int k = 0; k < coverage.cols(); ++k) {
+            double cov = coverage(f, k);
+            if (std::isnan(cov))
+                std::cout << "NaN ";
+            else
+                std::cout << cov << " ";
+        }
+        std::cout << "\n";
+    }
+
+    //Uniformity deviation analysis 
+    std::cout << "\n";
+    std::cout << "Facet        Max deviation (%)        Avg deviation (%)\n";
+    for (int f = 0; f < max_dev.size(); ++f) {
+        std::cout << std::setw(6) << f << " "
+                  << std::fixed << std::setprecision(2)
+                  << std::setw(18) << max_dev[f] << " "
+                  << std::setw(22) << avg_dev[f]
+                  << "\n";
+    }
+
+    return 0;
+}