diff --git a/benchmark/CMakeLists.txt b/benchmark/CMakeLists.txt new file mode 100644 index 000000000..85167e80a --- /dev/null +++ b/benchmark/CMakeLists.txt @@ -0,0 +1,81 @@ +cmake_minimum_required(VERSION 3.10) +project(VolestiBenchmark) + +set(CMAKE_CXX_STANDARD 17) +set(CMAKE_CXX_STANDARD_REQUIRED ON) + +if(NOT CMAKE_BUILD_TYPE) + set(CMAKE_BUILD_TYPE Release) +endif() +add_compile_options(-O3) + +add_compile_definitions(EIGEN_NO_DEBUG) +add_compile_definitions(BOOST_NO_AUTO_PTR) +add_compile_definitions(DISABLE_NLP_ORACLES) + +find_package(Boost REQUIRED COMPONENTS program_options) +find_package(Eigen3 3.3 REQUIRED NO_MODULE) + +find_library(LP_SOLVE NAMES liblpsolve55.so + PATHS /usr/lib/lp_solve /usr/local/lib) +find_path(LP_SOLVE_INCLUDE_DIR NAMES lp_lib.h + PATH_SUFFIXES lpsolve lp_solve) + +find_library(QD_LIB NAMES qd) + +if(NOT LP_SOLVE) + message(WARNING "lp_solve not found. ComputeInnerBall() might fail to link.") +else() + message(STATUS "Library lp_solve found: ${LP_SOLVE}") +endif() + +if(NOT LP_SOLVE_INCLUDE_DIR) + message(WARNING "lp_solve headers not found (lp_lib.h). Adjust include paths if needed.") +else() + message(STATUS "lp_solve headers found: ${LP_SOLVE_INCLUDE_DIR}") +endif() + +if(NOT QD_LIB) + message(WARNING "qd library not found. CRHMC / dd_real math will fail to link.") +else() + message(STATUS "Library qd found: ${QD_LIB}") +endif() + +include_directories( + ${CMAKE_CURRENT_SOURCE_DIR}/include + ${CMAKE_CURRENT_SOURCE_DIR}/../include + ${CMAKE_CURRENT_SOURCE_DIR}/../external + ${CMAKE_CURRENT_SOURCE_DIR}/../include/generators + ${CMAKE_CURRENT_SOURCE_DIR}/../include/convex_bodies + ${CMAKE_CURRENT_SOURCE_DIR}/../include/random_walks + ${CMAKE_CURRENT_SOURCE_DIR}/../include/preprocess +) + +if(LP_SOLVE_INCLUDE_DIR) + include_directories(${LP_SOLVE_INCLUDE_DIR}) +endif() + +set(SOURCES + src/benchmark_main.cpp + src/benchmark_run.cpp + src/walk_parameters.cpp + src/benchmark_utils.cpp + src/walk_result.cpp + src/walk_registry.cpp +) + +add_executable(benchmark_run ${SOURCES}) + +target_link_libraries(benchmark_run + Boost::program_options + Eigen3::Eigen + m +) + +if(LP_SOLVE) + target_link_libraries(benchmark_run ${LP_SOLVE}) +endif() + +if(QD_LIB) + target_link_libraries(benchmark_run ${QD_LIB}) +endif() \ No newline at end of file diff --git a/benchmark/Readme.md b/benchmark/Readme.md new file mode 100644 index 000000000..63050e341 --- /dev/null +++ b/benchmark/Readme.md @@ -0,0 +1,100 @@ +# Volesti Benchmark Suite + +This project is a benchmarking suite for testing various random walk algorithms on high-dimensional convex polytopes found inside the volesti library + +## Usage + +Create a build folder, and run Cmake and Make. Make sure all the needed libraries are installed. + +```bash +mkdir build && cd build +cmake .. +make +``` + +Then, just run the executable: + +```bash +./benchmark_run +``` + +## Command-Line Arguments + +While not necessary, the executable accepts the following optional parameters + +### `-h` or `--help` + +Prints the help message and available options, then exits safely. + +--- + +### `-c ` or `--config ` + +Specify a custom path to your JSON configuration file. + +- **Default:** `../config/walk_config.json` + +--- + +### `-d ` or `--dim ` + +Set the dimension for generating standard polytopes (e.g., Cube, Simplex). + +- **Default:** The `"dimension"` value in your JSON config. +- **Note:** *Custom Polytopes* will use its own dimesnion based on the csv files. + +--- + +### `-p ` or `--polytope ` + +Choose the mathematical shape of the polytope. + +- **Valid Options:** Cube, Simplex, Birkhoff, Cross, OrderPolytope, Custom + +--- + +### `-w ` or `--walk ` + +Specify exactly which random walk to run. Walk names are case-sensitive. + +- **Valid Options:** + All, BallWalk, BilliardWalk, AcceleratedBilliardWalk, SparseBilliardWalk, CDHRWalk, RDHRWalk, DikinWalk, JohnWalk, VaidyaWalk, GaussianBallWalk, GaussianCDHRWalk, BilliardShakeAndBakeWalk, ShakeAndBakeWalk, BCDHRWalk, BRDHRWalk + | + +## Examples + +### Run the default benchmark + +(Uses the configuration file to determine the experiment parameters) + +```bash +./benchmark_run +``` + +### Test a specific algorithm at a specific dimension and a specific Polytope + +(Runs only the Accelerated Billiard Walk in 25 dimensions for the Simplex polytope) + +```bash +./benchmark_run -p Simplex -w AcceleratedBilliardWalk -d 25 +``` + +### Use your custom csv polytope + +```bash +./benchmark_run -p Custom -w BilliardWalk +``` + +### 3. Use shorthand flags for a custom configuration + +(Runs the Billiard Walk in 100 dimensions using a custom JSON file) + +```bash +./benchmark_run -w BilliardWalk -d 100 -c ../config/walk_config.json +``` + +You actually don't have to use any argument at all. Just run the "benchmark_run" file and it will retrieve all the necessary parameters from the configuration file. + +## Configuration (JSON) + +Algorithm-specific parameters are managed entirely via the JSON configuration file. To change them, simply edit `config/walk_config.json`. diff --git a/benchmark/config/available_choices.md b/benchmark/config/available_choices.md new file mode 100644 index 000000000..4cf5bff57 --- /dev/null +++ b/benchmark/config/available_choices.md @@ -0,0 +1,76 @@ +# Available Choices + +## Target ESS + +The code will use a dynamic batch size in order to sample the necessary number of samples needed to reach the targert ESS + +## Time limit + +The maximum time the code runs. If the time runs out, the code will stop at the first available opportunity. + +## Base seed + +The seed used by the random number generator. Use any number you like. + +## Dimension + +The number of dimesnion the generated polytope will have. This only applies to polytopes generated by the generator functions. It is ignored for custom polytopes. + +## Rotation angle + +You can rotate the input polytope by this number over all pairs of dimensions. Angle is in degrees. Use 0 to avoid rotation. + +## Polytope Choice + +Choose the polytope that will be sampled from. +Available polytope choices are: Cube, Simplex, Birkhoff, Cross, OrderPolytope, Custom + +## Custom A file and Custom b file + +These are the names of the A and b files of your custom polytope defining the polytope Ax<=b. Csv files are expected, with no labels. +Place them inside build folder next to the executable. + +## Dynamic batch size + +If you have this choice to true, the code will run a dynamic batch size in order to hit target ESS. +If you have this choice to false, the code will generate exactly "samples" amount of samples, the parameter defined for each walk method. + +## Write to file + +Set this to true in order to create txt files holding all the generated samples. Useful for later use of the samples. + +## Rounding + +If this is set to true, the code will first call a rounding function on the input polytope. +The sampling will take place on the rounded polytope before the sample are reverted back to the original and returned. + +## Rounding_method + +Selects the rounding method. Available choices are: + +- max_ellipsoid +- log_barrier +- vaidya_barrier +- volumetric_barrier + +## Auto walk + +If set to true, the code will decide what algorithm to use based on the input polytope. +Right now it decides only based on dimension, but will be updated in the future. + +## Show console logs + +Keep your console clean by removing most of the cosnole messages. You will still see the progress bar and final a few mesages. + +## Show menu + +If set to true, a menu will show up when you run the code. The menu will ask you to choose all the parameters you want for your experiment or show extended help. Useful for new users. + +## Walks + +For each walk, you can change these parameters: + +- enabled: This will determine if this method will run or not when sampling +- samples: How many samples the code will take if dynamic batch size is also false. If dynamic batch size is on, target ESS will be used instead to determine when the execution stops. +- walk_len: It is a thinning parameter. The code will save the sample only after walk_len samples are sampled. +walk_len = walk_base + dim * walk_multiplier" diff --git a/benchmark/config/walk_config.json b/benchmark/config/walk_config.json new file mode 100644 index 000000000..05392722e --- /dev/null +++ b/benchmark/config/walk_config.json @@ -0,0 +1,37 @@ +{ + "global_settings": { + "target_ESS": 3000, + "time_limit_sec": 30.0, + "base_seed": 42, + "dimensions": [50], + "rotation_angle": 53, + "polytope_choice": "Cube", + "custom_A_file": "Cube_50_A.csv", + "custom_b_file": "Cube_50_b.csv", + "dynamic_batch_size": true, + "write_to_file": true, + "rounding": true, + "rounding_method": "max_ellipsoid", + "auto_walk": false, + "show_console_logs": true, + "show_menu": true + }, +"walks": { + "BallWalk": { "enabled": true, "samples": 20000, "walk_len_multiplier": 2, "walk_len_base": 0 }, + "BilliardWalk": { "enabled": false, "samples": 2000, "walk_len_multiplier": 0, "walk_len_base": 1 }, + "AcceleratedBilliardWalk": { "enabled": false, "samples": 2000, "walk_len_multiplier": 0, "walk_len_base": 1 }, + "SparseBilliardWalk": { "enabled": false, "samples": 1000, "walk_len_multiplier": 0, "walk_len_base": 1 }, + "CDHRWalk": { "enabled": false, "samples": 1000, "walk_len_multiplier": 2, "walk_len_base": 0 }, + "RDHRWalk": { "enabled": false, "samples": 1000, "walk_len_multiplier": 2, "walk_len_base": 0 }, + "DikinWalk": { "enabled": false, "samples": 1000, "walk_len_multiplier": 2, "walk_len_base": 0 }, + "JohnWalk": { "enabled": false, "samples": 1000, "walk_len_multiplier": 2, "walk_len_base": 0 }, + "VaidyaWalk": { "enabled": false, "samples": 1000, "walk_len_multiplier": 2, "walk_len_base": 0 }, + "GaussianBallWalk": { "enabled": false, "samples": 1000, "walk_len_multiplier": 2, "walk_len_base": 0, "a_i_param": 1.0 }, + "GaussianCDHRWalk": { "enabled": false, "samples": 1000, "walk_len_multiplier": 2, "walk_len_base": 0, "a_i_param": 1.0 }, + "BilliardShakeAndBakeWalk": { "enabled": false, "samples": 1000, "walk_len_multiplier": 2, "walk_len_base": 0 }, + "ShakeAndBakeWalk": { "enabled": false, "samples": 1000, "walk_len_multiplier": 2, "walk_len_base": 0 }, + "BCDHRWalk": { "enabled": false, "samples": 1000, "walk_len_multiplier": 2, "walk_len_base": 0 }, + "BRDHRWalk": { "enabled": false, "samples": 1000, "walk_len_multiplier": 2, "walk_len_base": 0 }, + "CRHMCWalk": { "enabled": false, "samples": 1000, "walk_len_multiplier": 0, "walk_len_base": 1 } + } +} \ No newline at end of file diff --git a/benchmark/include/benchmark_run.hpp b/benchmark/include/benchmark_run.hpp new file mode 100644 index 000000000..49715fb08 --- /dev/null +++ b/benchmark/include/benchmark_run.hpp @@ -0,0 +1,3 @@ +#pragma once + +int run_benchmark(int argc, char** argv); \ No newline at end of file diff --git a/benchmark/include/benchmark_utils.hpp b/benchmark/include/benchmark_utils.hpp new file mode 100644 index 000000000..10f2ca563 --- /dev/null +++ b/benchmark/include/benchmark_utils.hpp @@ -0,0 +1,46 @@ +#pragma once + +#include "core_types.hpp" +#include +#include +#include + +//In this file you can find the definition of supportive functions. +//Timer is a class to count the time for each random walk. +//write to file is used to write the results to a file. +//vector to eigen converts a vector in an eigen compatible form. (eigen is numpy for c++) +//determine auto walk is used to choose a walk based on the polytope structure. + + +class Timer { +public: + Timer(const std::string& name = ""); + void start(); + double stop(const std::string& label = ""); + double get_total_time() const; + +private: + std::string walk_name; + std::chrono::steady_clock::time_point start_time; + double total_time; + bool is_running; +}; + + +void write_to_file(std::string filename, std::vector const& randPoints); + + +template +MT vector_to_eigen(const std::vector& someSamples) { + if (someSamples.empty()) { + return MT(); // Return empty matrix if no samples + } + + MT samples(someSamples[0].dimension(), someSamples.size()); + for (size_t jj = 0; jj < someSamples.size(); ++jj) { + samples.col(jj) = someSamples[jj].getCoefficients(); + } + return samples; +} + +std::string determine_auto_walk(unsigned int dim); \ No newline at end of file diff --git a/benchmark/include/core_types.hpp b/benchmark/include/core_types.hpp new file mode 100644 index 000000000..b4366c3aa --- /dev/null +++ b/benchmark/include/core_types.hpp @@ -0,0 +1,33 @@ +#pragma once + +#include +#include +#include + +#include +#include "Eigen/Eigen" + +#include "cartesian_geom/cartesian_kernel.h" +#include "convex_bodies/hpolytope.h" +#include "random_walks/random_walks.hpp" + +// In this file you can find the core type definitions and commonly used aliases +// for the project. +// +// It defines: +// Numeric type (NT) used across computations. +// Geometric kernel (Cartesian) and Point representation. +// Matrix (MT) and vector (VT) types using Eigen (similar to NumPy in C++). +// Random number generator type based on Boost. +// HPolytope type used to represent convex bodies. + +typedef double NT; +typedef Cartesian Kernel; +typedef Kernel::Point Point; +typedef Eigen::Matrix MT; +typedef Eigen::Matrix VT; +typedef BoostRandomNumberGenerator RNGType; +typedef HPolytope HPOLYTOPE; + +//it is only initialized once inside src/benchmark_utils.cpp +extern PushBackWalkPolicy push_back_policy; \ No newline at end of file diff --git a/benchmark/include/diagnostics.hpp b/benchmark/include/diagnostics.hpp new file mode 100644 index 000000000..b8a7d5065 --- /dev/null +++ b/benchmark/include/diagnostics.hpp @@ -0,0 +1,64 @@ +#pragma once + +#include +#include +#include + +#include "core_types.hpp" +#include "benchmark_utils.hpp" + +#include "sampling/sample_correlation_matrices.hpp" +#include "matrix_operations/EigenvaluesProblems.h" +#include "diagnostics/effective_sample_size.hpp" +#include "diagnostics/univariate_psrf.hpp" +#include "diagnostics/scaling_ratio.hpp" +#include "diagnostics/KS_test.hpp" + +// In this file you can find utility functions for computing diagnostic +// metrics on samples produced by random walks. + +// The file works with Eigen matrices and uses helper utilities such as +// vector_to_eigen to convert sample containers into matrix form. + +// Computes ESS on a given Eigen matrix and returns the min ESS +template +unsigned int compute_ess(const MT& samples) { + unsigned int min_ess = 0; + // call internal ESS function + VT ess_vector = effective_sample_size(samples, min_ess); + return min_ess; +} + +// Computes PSRF for all accumulated points +template +double compute_psrf(const std::vector& someSamples) { + // Convert the vector of points to an Eigen matrix using our utility + MT finalSamples = vector_to_eigen(someSamples); + + // Call Volesti's internal PSRF function + VT psrf = univariate_psrf(finalSamples); + return psrf.maxCoeff(); +} + + +// Struct to neatly pass the KS results back +struct KSTestResult { + double ks_stat; + double p_val; + std::vector observed; + std::vector expected; +}; + +// Computes the KS Test and scaling ratios +template +KSTestResult compute_ks_test(const PolytopeType& polytope, const MT& samples_mat, double current_ESS) { + // Calculate a safe thinning factor based on ESS + double safe_ess = (current_ESS > 0) ? current_ESS : 1.0; + int computed_thin = static_cast(samples_mat.cols() / safe_ess); + int thin_factor = std::max(10, computed_thin * 2); + + // Call KS test + auto [ks_stat, p_val, observed, expected] = global_scaling_test(polytope, samples_mat, thin_factor); + + return {ks_stat, p_val, observed, expected}; +} \ No newline at end of file diff --git a/benchmark/include/dynamic_batch_size.hpp b/benchmark/include/dynamic_batch_size.hpp new file mode 100644 index 000000000..5285ec863 --- /dev/null +++ b/benchmark/include/dynamic_batch_size.hpp @@ -0,0 +1,114 @@ +#pragma once + +#include + +// Checks if a number is "7-smooth" (only divisible by 2, 3, 5, or 7) +inline bool is_fft_friendly(unsigned int n) { + if (n == 0) return false; + while (n % 2 == 0) n /= 2; + while (n % 3 == 0) n /= 3; + while (n % 5 == 0) n /= 5; + while (n % 7 == 0) n /= 7; + return n == 1; // If it reduces exactly to 1, the FFT will process it instantly. +} + +// Finds the nearest FFT-friendly number by bumping it up by 1 until it fits +inline unsigned int get_next_fft_friendly_size(unsigned int target) { + unsigned int n = target; + while (!is_fft_friendly(n)) { + n++; + } + return n; +} + +// Calculates the optimal next batch size based on the current ESS +inline unsigned int compute_next_batch_size( + unsigned int target_ESS, + unsigned int current_ESS, + unsigned int previous_ESS, + size_t total_samples, + unsigned int current_batch_size, + unsigned int dimension, + double remaining_time_sec, + double samples_per_sec) +{ + unsigned int base_cap = std::max(500000u, target_ESS * 10); + unsigned int penalty = 1000 * dimension; + unsigned int MAX_BATCH_SIZE = (base_cap > penalty) ? (base_cap - penalty) : 1000; + MAX_BATCH_SIZE = std::max(1000u, std::min(MAX_BATCH_SIZE, 2000000u)); + + double ess_per_sample = static_cast(current_ESS) / static_cast(total_samples); + + // The next lines attempt to catch a stuck sampler by checking the ESS efficiency between batches + // calculate how efficiently this specific batch generated ESS + unsigned int ess_gained = (current_ESS > previous_ESS) ? (current_ESS - previous_ESS) : 0; + double marginal_ess_per_sample = static_cast(ess_gained) / static_cast(current_batch_size); + size_t previous_samples = (total_samples > current_batch_size) ? (total_samples - current_batch_size) : 0; + double previous_ess_per_sample = (previous_samples > 0) ? (static_cast(previous_ESS) / static_cast(previous_samples)) : 0.0; + + // warning check + if (previous_samples > 0) { + // Condition 1: efficiency dropped by 90% or more compared to historical average + if (previous_ess_per_sample > 1e-6 && marginal_ess_per_sample < (previous_ess_per_sample * 0.1)) { + std::cerr << " | [WARNING] Sampler might be stuck! ESS yield dropped heavily." << std::flush; + } + // Condition 2: absolute terrible mixing (>10,000 samples per 1 ESS) + else if (marginal_ess_per_sample < 1e-4) { + std::cerr << " | [WARNING] Terrible mixing detected. Sampler may be stuck in a corner." << std::flush; + } + } + + unsigned int next_batch_size = 0; + unsigned int min_viable_batch = 2500u; + + // The main idea is to ask for samples based on how many samples we need for 1 ESS. + if (ess_per_sample > 1e-6) { + unsigned int remaining_ESS = target_ESS - current_ESS; + + // we always ask for enough points to generate at least ~15 ESS. + min_viable_batch = static_cast(15.0 / ess_per_sample); + + // Keep an absolute minimum, but never let the minimum exceed the calculated maximum + min_viable_batch = std::max(min_viable_batch, 2500u); + min_viable_batch = std::min(min_viable_batch, MAX_BATCH_SIZE); + + if (current_ESS < (target_ESS * 0.85)) { + // we are less than 85% to the target. + unsigned int samples_needed = static_cast((remaining_ESS / ess_per_sample) * 0.80); + next_batch_size = std::max(samples_needed, min_viable_batch); + } else { + // we are in the final 15%. + unsigned int samples_needed = static_cast((remaining_ESS / ess_per_sample) * 1.10); + next_batch_size = std::max(samples_needed, min_viable_batch); + } + } + else { + // Fallback if ess_per_sample is virtually zero (terrible mixing) + // Back off the aggressive multiplier. + next_batch_size = static_cast(current_batch_size * 1.25); + + // If we are stuck, do not let it grow all the way to MAX_BATCH_SIZE. + // Cap the runaway growth at 25% of our maximum (with a floor of 1000). + unsigned int safe_stuck_cap = std::max(1000u, MAX_BATCH_SIZE / 4); + next_batch_size = std::min(next_batch_size, safe_stuck_cap); + } + + if (remaining_time_sec > 0 && samples_per_sec > 0) { + unsigned int time_budget_batch = static_cast(remaining_time_sec * samples_per_sec * 1.2); + next_batch_size = std::min(next_batch_size, std::max(time_budget_batch, min_viable_batch)); + } + + // Calculate the raw batch size bounded by our maximums + unsigned int raw_next_batch = std::min(next_batch_size, MAX_BATCH_SIZE); + + if (raw_next_batch == 0) return 0; + + // We calculate what the absolute total number of samples will be after this batch + unsigned int target_total_samples = total_samples + raw_next_batch; + + // We bump the total up slightly (usually < 100 points) to the next Smooth Number + unsigned int optimal_total_samples = get_next_fft_friendly_size(target_total_samples); + + // Return the adjusted batch size + return optimal_total_samples - total_samples; +} \ No newline at end of file diff --git a/benchmark/include/geometry_utils.hpp b/benchmark/include/geometry_utils.hpp new file mode 100644 index 000000000..56dc79c0f --- /dev/null +++ b/benchmark/include/geometry_utils.hpp @@ -0,0 +1,41 @@ +#pragma once + +#include "core_types.hpp" +#include + +// Useful function to rotate polytope +template +HPOLYTOPE rotate_all_dims(const HPOLYTOPE& P, typename HPOLYTOPE::NT angle) +{ + using NT = typename HPOLYTOPE::NT; + int dim = P.dimension(); + + // Build global rotation matrix + Eigen::Matrix R = + Eigen::Matrix::Identity(dim, dim); + + NT c = std::cos(angle); + NT s = std::sin(angle); + + // Apply rotation in each adjacent coordinate plane + for (int k = 0; k < dim - 1; ++k) { + Eigen::Matrix Rk = + Eigen::Matrix::Identity(dim, dim); + + Rk(k, k) = c; + Rk(k, k+1) = -s; + Rk(k+1, k) = s; + Rk(k+1, k+1) = c; + + R = R * Rk; // compose rotations + } + + // Extract A and b + auto A = P.get_mat(); + auto b = P.get_vec(); + + // Apply A' = A R + Eigen::Matrix A_rot = A * R; + + return HPOLYTOPE(dim, A_rot, b); +} \ No newline at end of file diff --git a/benchmark/include/menu.hpp b/benchmark/include/menu.hpp new file mode 100644 index 000000000..90ed9dd3b --- /dev/null +++ b/benchmark/include/menu.hpp @@ -0,0 +1,489 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include "walk_parameters.hpp" + +// Helper function to safely get an integer or quit via 'q'/'Q' from ANY prompt +inline int get_valid_int(const std::string& prompt, int min_val, int max_val) { + std::string input; + while (true) { + std::cout << prompt; + if (!(std::cin >> input)) { + std::cin.clear(); + std::cin.ignore(std::numeric_limits::max(), '\n'); + std::cout << "[!] Invalid input. Enter a number or 'q' to quit: "; + continue; + } + + // Check for exit command anywhere + if (input == "q" || input == "Q") { + std::cout << "\nExiting program. Goodbye!\n"; + std::exit(0); + } + + // Try parsing the string into an integer safely + try { + size_t idx; + int choice = std::stoi(input, &idx); + + // Ensure the entire string was consumed (rejects strings like "5abc") + if (idx != input.size()) { + std::cout << "[!] Invalid format. Please enter a valid number or 'q' to quit.\n"; + continue; + } + + if (choice < min_val || choice > max_val) { + std::cout << "[!] Out of range. Please choose between " << min_val << " and " << max_val << " (or 'q' to quit).\n"; + continue; + } + return choice; + } catch (const std::exception&) { + std::cout << "[!] Invalid input. Please enter a valid number or 'q' to quit.\n"; + } + } +} + +// Helper function to get a comma-separated list of integers +inline std::vector get_valid_dimension_list(const std::string& prompt) { + std::string input; + while (true) { + std::cout << prompt; + + // std::ws skips any leftover newlines in the input buffer from previous prompts + if (!std::getline(std::cin >> std::ws, input)) { + std::cin.clear(); + std::cout << "[!] Invalid input.\n"; + continue; + } + + if (input == "q" || input == "Q") { + std::cout << "\nExiting program. Goodbye!\n"; + std::exit(0); + } + + std::vector dims; + std::stringstream ss(input); + std::string token; + bool error = false; + + // Split the string by commas + while (std::getline(ss, token, ',')) { + // Strip any accidental spaces (e.g. if they type "10, 20, 30") + token.erase(std::remove_if(token.begin(), token.end(), ::isspace), token.end()); + + if (token.empty()) continue; + + try { + size_t idx; + int val = std::stoi(token, &idx); + + // Ensure it's a valid integer within a reasonable range + if (idx != token.size() || val < 1 || val > 100000) { + error = true; + break; + } + dims.push_back(static_cast(val)); + } catch (...) { + error = true; + break; + } + } + + if (error || dims.empty()) { + std::cout << "[!] Invalid format. Enter positive integers separated by commas (e.g., 10,20,30).\n"; + continue; + } + + return dims; + } +} + +inline void show_help() { + int choice = -1; + + while (true) { + // Clear the console screen (cross-platform compatible) +#if defined(_WIN32) + int ret = std::system("cls"); +#else + int ret = std::system("clear"); +#endif + (void)ret; // Silence the "ignoring return value" warning + + // Print the fixed menu at the top + std::cout << "\n" << std::string(60, '-') << "\n"; + std::cout << "*************************** HELP ***************************\n"; + std::cout << std::string(60, '-') << "\n"; + std::cout << "Select a topic to learn more about:\n"; + std::cout << " 1. What is an MCMC method?\n"; + std::cout << " 2. What is ESS (Effective Sample Size)?\n"; + std::cout << " 3. What is PSRF?\n"; + std::cout << " 4. What is the KS test?\n"; + std::cout << " 5. Polytope Rounding\n"; + std::cout << " 6. Geometric Methods (Ball, Billiards, Hit-and-Run)\n"; + std::cout << " 7. Barrier Methods (Dikin, John, Vaidya, CRHMC)\n"; + std::cout << " 8. Boundary & Other Methods (Shake & Bake)\n"; + std::cout << " 0. Return to Main Menu\n"; + + std::cout << "\n" << std::string(60, '=') << "\n"; + + // Display the selected information block + if (choice == -1) { + std::cout << "Please select an option from the menu above.\n"; + } else { + switch (choice) { + case 1: + std::cout << R"(--- Markov Chain Monte Carlo (MCMC) --- +MCMC is a class of algorithms used to sample from a probability +distribution. Because high-dimensional polytopes are difficult +to sample from directly, MCMC constructs a random walk (a +Markov chain) where each step depends only on the previous one. +Over time, the steps of this walk distribute themselves +according to the target distribution (usually uniform for +polytopes).)" << '\n'; + break; + case 2: + std::cout << R"(--- Effective Sample Size (ESS) --- +Since MCMC samples are generated via a random walk, consecutive +samples are correlated (they are near each other in space). +ESS calculates how many statistically independent samples your +autocorrelated chain is equivalent to. + +For example, 10,000 raw samples might only yield an ESS of 500, +meaning the walk is moving slowly through the space. A higher +ESS indicates better mixing and a more efficient algorithm.)" << '\n'; + break; + case 3: + std::cout << R"(--- Potential Scale Reduction Factor (PSRF) --- +Also known as the Gelman-Rubin diagnostic. It is used to evaluate +if your MCMC walk has converged to the target distribution. + +It compares the variance between multiple independent random +walks to the variance within those individual walks. A PSRF +value close to 1.0 (typically < 1.1) suggests that the chains +have converged and are exploring the same space properly.)" << '\n'; + break; + case 4: + std::cout << R"(--- Kolmogorov-Smirnov (KS) Test --- +The KS test is a statistical test used to compare a sample +distribution with a reference probability distribution, or to +compare two sample distributions. + +In this benchmark, it is often used to check the 1-dimensional +marginals of the sampled points to see if they match the +theoretically expected uniform distribution of the polytope, +verifying the correctness of the walk.)" << '\n'; + break; + case 5: + std::cout << R"(--- Polytope Rounding --- +Polytopes can be highly skewed, elongated, or "skinny". Random +walks struggle in skinny spaces because they easily get stuck +bouncing off the walls, leading to high autocorrelation. + +Rounding computes an ellipsoid (e.g., Maximum Volume Inscribed +Ellipsoid) that fits inside the polytope, and uses it to apply +a linear transformation. This transformation makes the polytope +look more like a sphere (isotropic). The random walk runs much +faster in this rounded space, and the generated samples are then +transformed back to the original space.)" << '\n'; + break; + case 6: + std::cout << R"(--- Geometric Methods --- +These methods rely purely on the geometric properties of the space: +- Hit-and-Run: Picks a random direction, computes the line segment + inside the polytope along that direction, and jumps to a uniform + random point on that segment. +- Ball Walk: Proposes a point uniformly within a small ball around + the current point. Rejects the step if it falls outside the polytope. +- Billiards: Simulates particle trajectories bouncing off the + polytope's inner boundaries, like a billiard ball.)" << '\n'; + break; + case 7: + std::cout << R"(--- Barrier Methods --- +These methods utilize barrier functions (often borrowed from +Interior Point methods in optimization) to keep the walk away from +the boundaries: +- Dikin Walk: Uses the log-barrier Hessian to define an ellipsoid + around the current point that is guaranteed to be inside the + polytope, proposing the next step from within this ellipsoid. +- Vaidya / John Walks: Use more complex barrier functions + (Vaidya's barrier or John's ellipsoid) to take larger steps + without hitting the boundary. +- CRHMC: Simulates a particle moving inside the polytope under + the Hamiltonian laws of motion, where the boundaries act as + immovable walls.)" << '\n'; + break; + case 8: + std::cout << R"(--- Boundary & Other Methods --- +While standard methods sample the interior volume of the polytope, +some tasks require sampling exactly from the surface (the facets). + +- Shake & Bake: An algorithm specifically designed to generate + points uniformly distributed on the boundary of a polytope. It + "shakes" to find a direction and "bakes" by moving along the + boundary facets.)" << '\n'; + break; + } + } + std::cout << std::string(60, '=') << "\n"; + + // Prompt the user for the next action + choice = get_valid_int("Enter your choice (0-8): ", 0, 8); + + if (choice == 0) { + // Clear the screen one last time before returning to the main menu +#if defined(_WIN32) + int final_ret = std::system("cls"); +#else + int final_ret = std::system("clear"); +#endif + (void)final_ret; // Silence the warning + return; + } + } +} + +// Handles the setup flow when the user chooses Option 1 +inline void setup_benchmark_options(BenchmarkConfig& config, std::string& walk_choice) { + std::cout << "\n--- Benchmark Setup (Type 'q' to quit at any prompt) ---\n"; + + // 1. Polytope Choice + std::string poly_prompt = "Select a Polytope:\n" + "1. Cube\n" + "2. Simplex\n" + "3. Birkhoff\n" + "4. Cross\n" + "5. OrderPolytope\n" + "6. Custom\n" + "Choice (1-6): "; + + int p_choice = get_valid_int(poly_prompt, 1, 6); + switch (p_choice) { + case 1: config.polytope_choice = "Cube"; break; + case 2: config.polytope_choice = "Simplex"; break; + case 3: config.polytope_choice = "Birkhoff"; break; + case 4: config.polytope_choice = "Cross"; break; + case 5: config.polytope_choice = "OrderPolytope"; break; + case 6: + config.polytope_choice = "Custom"; + std::cout << "\n[Note] You selected Custom. The tool will use the 'custom_A_file'\n" + << "and 'custom_b_file' exactly as defined in your JSON config.\n"; + break; + } + + // 2. Dimension (Skip if Custom Polytope) + if (config.polytope_choice != "Custom") { + std::cout << "\n"; + // Directly overwrite the JSON dimensions array with the user's list + config.dimensions = get_valid_dimension_list("Enter dimensions separated by commas (e.g., 10,20,50): "); + } else { + config.dimensions.clear(); + config.dimensions.push_back(0); + } + + // 3. Stopping Criterion + std::cout << "\n"; + std::string stop_prompt = "Select stopping criterion:\n" + "1. Target ESS\n" + "2. Fixed number of samples\n" + "Choice (1-2): "; + int stop_choice = get_valid_int(stop_prompt, 1, 2); + + if (stop_choice == 1) { + config.use_dynamic_batch = true; + config.target_ESS = get_valid_int("\nEnter Target ESS (10 to 100000): ", 10, 100000); + } else { + config.use_dynamic_batch = false; + int fixed_samples = get_valid_int("\nEnter number of samples to generate (1 to 100000): ", 1, 100000); + + // Override the 'samples' value for all walks in the configuration. + // This ensures your get_initial_batch_size() function reads the correct user input. + for (auto& pair : config.walk_settings) { + pair.second.samples = fixed_samples; + } + } + + // 4. Time Limit + std::cout << "\n"; + int time_lim = get_valid_int("Enter time limit in seconds (1 to 86400): ", 1, 86400); + config.time_limit_sec = static_cast(time_lim); + + // 5. Base Seed + std::cout << "\n"; + config.base_seed = get_valid_int("Enter base seed (e.g., 42): ", 0, 2147483647); + + // 6. Rotation Angle + std::cout << "\n"; + config.angle = get_valid_int("Enter rotation angle in degrees (0 to 360): ", 0, 360); + + // 7. Write to File + std::cout << "\n"; + int write_choice = get_valid_int("Write samples to file?\n1. Yes\n2. No\nChoice (1-2): ", 1, 2); + config.write_to_file = (write_choice == 1); + + // 8. Rounding + std::cout << "\n"; + std::cout << "Select a rounding method:\n"; + std::cout << "1. None (Disabled)\n"; + std::cout << "2. Max Ellipsoid\n"; + std::cout << "3. Log Barrier\n"; + std::cout << "4. Vaidya Barrier\n"; + std::cout << "5. Volumetric Barrier\n"; + int r_choice = get_valid_int("Choice (1-5): ", 1, 5); + + switch (r_choice) { + case 1: + config.rounding = false; + config.rounding_method = "none"; + break; + case 2: + config.rounding = true; + config.rounding_method = "max_ellipsoid"; + break; + case 3: + config.rounding = true; + config.rounding_method = "log_barrier"; + break; + case 4: + config.rounding = true; + config.rounding_method = "vaidya_barrier"; + break; + case 5: + config.rounding = true; + config.rounding_method = "volumetric_barrier"; + break; + } + // 9. Method Choice + std::cout << "\n"; + std::string strategy_prompt = "Select a Walk Method Strategy:\n" + "1. All (Run all enabled in JSON)\n" + "2. Auto (Select based on dimension)\n" + "3. Manual Selection\n" + "Choice (1-3): "; + + int strat_choice = get_valid_int(strategy_prompt, 1, 3); + + if (strat_choice == 1) { + walk_choice = "All"; + config.auto_walk = false; + } + else if (strat_choice == 2) { + walk_choice = "All"; + config.auto_walk = true; + } + else { + config.auto_walk = false; + std::string dist_prompt = "\nSelect Target Distribution:\n" + "1. Uniform\n" + "2. Exponential (Gaussian)\n" + "Choice (1-2): "; + int dist_choice = get_valid_int(dist_prompt, 1, 2); + + if (dist_choice == 2) { + std::string exp_prompt = "\nSelect Exponential (Gaussian) Method:\n" + "1. Gaussian Ball\n" + "2. Gaussian Coordinate Direction Hit and Run\n" + "Choice (1-2): "; + int exp_choice = get_valid_int(exp_prompt, 1, 2); + switch(exp_choice) { + case 1: walk_choice = "GaussianBallWalk"; break; + case 2: walk_choice = "GaussianCDHRWalk"; break; + } + } + else { + std::string cat_prompt = "\nSelect Uniform Method Category:\n" + "1. Geometric (Ball, Billiards, Hit-and-Run)\n" + "2. Barrier (Dikin, John, Vaidya, CRHMC)\n" + "3. Boundary & Other (Shake & Bake)\n" + "Choice (1-3): "; + int cat_choice = get_valid_int(cat_prompt, 1, 3); + + if (cat_choice == 1) { + std::string geom_prompt = "\nSelect Geometric Method:\n" + "1. Ball\n" + "2. Billiard\n" + "3. Accelerated Billiard\n" + "4. Sparse Billiard\n" + "5. Coordinate Direction Hit and Run (CDHR)\n" + "6. Random Direction Hit and Run (RDHR)\n" + "Choice (1-6): "; + int geom_choice = get_valid_int(geom_prompt, 1, 6); + switch(geom_choice) { + case 1: walk_choice = "BallWalk"; break; + case 2: walk_choice = "BilliardWalk"; break; + case 3: walk_choice = "AcceleratedBilliardWalk"; break; + case 4: walk_choice = "SparseBilliardWalk"; break; + case 5: walk_choice = "CDHRWalk"; break; + case 6: walk_choice = "RDHRWalk"; break; + } + } + else if (cat_choice == 2) { + std::string bar_prompt = "\nSelect Barrier Method:\n" + "1. Dikin\n" + "2. John\n" + "3. Vaidya\n" + "4. Constrained Riemannian Hamiltonian Monte Carlo (CRHMC)\n" + "Choice (1-4): "; + int bar_choice = get_valid_int(bar_prompt, 1, 4); + switch(bar_choice) { + case 1: walk_choice = "DikinWalk"; break; + case 2: walk_choice = "JohnWalk"; break; + case 3: walk_choice = "VaidyaWalk"; break; + case 4: walk_choice = "CRHMCWalk"; break; + } + } + else { + std::string other_prompt = "\nSelect Boundary & Other Method:\n" + "1. Billiard Shake and Bake\n" + "2. Shake and Bake\n" + "3. Boundary Coordinate Direction Hit and Run\n" + "4. Boundary Random Direction Hit and Run\n" + "Choice (1-4): "; + int other_choice = get_valid_int(other_prompt, 1, 4); + switch(other_choice) { + case 1: walk_choice = "BilliardShakeAndBakeWalk"; break; + case 2: walk_choice = "ShakeAndBakeWalk"; break; + case 3: walk_choice = "BCDHRWalk"; break; + case 4: walk_choice = "BRDHRWalk"; break; + } + } + } + } + + std::cout << "\nSetup complete! Proceeding to benchmark...\n"; +} + +// Returns true if the benchmark should proceed, false if the user chose to quit. +inline bool run_interactive_menu(BenchmarkConfig& config, std::string& walk_choice) { + std::cout << "\n==========================================\n"; + std::cout << " Welcome to the Volesti Benchmark Tool \n"; + std::cout << "==========================================\n\n"; + + while (true) { + std::cout << "Main Menu:\n"; + std::cout << "1. Start benchmark\n"; + std::cout << "2. Show help\n"; + std::cout << "3. Quit\n"; + + int choice = get_valid_int("Enter your choice (1-3): ", 1, 3); + + switch (choice) { + case 1: + setup_benchmark_options(config, walk_choice); + return true; + case 2: + show_help(); + break; + case 3: + std::cout << "\nExiting program. Goodbye!\n"; + return false; + } + } +} \ No newline at end of file diff --git a/benchmark/include/polytope_generation.hpp b/benchmark/include/polytope_generation.hpp new file mode 100644 index 000000000..bc2bb8437 --- /dev/null +++ b/benchmark/include/polytope_generation.hpp @@ -0,0 +1,69 @@ +#pragma once + +#include "core_types.hpp" +#include "walk_parameters.hpp" + +#include +#include +#include +#include +#include +#include + +#include "known_polytope_generators.h" +#include "order_polytope_generator.h" +#include "custom_generators.h" + +using PolytopeGeneratorFn = std::function; + +inline HPOLYTOPE create_polytope(const std::string& choice, unsigned int dim, const BenchmarkConfig& config) { + + static std::map factory = { + {"Cube", [](unsigned int d, const BenchmarkConfig&) { + return generate_cube(d, false); + }}, + {"Simplex", [](unsigned int d, const BenchmarkConfig&) { + return generate_simplex(d, false); + }}, + {"Birkhoff", [](unsigned int d, const BenchmarkConfig&) { + return generate_birkhoff(d); + }}, + {"Cross", [](unsigned int d, const BenchmarkConfig&) { + return generate_cross(d, false); + }}, + {"OrderPolytope", [](unsigned int d, const BenchmarkConfig& cfg) { + unsigned int m = 3 * d; + int seed = cfg.base_seed + d; + return random_orderpoly(d, m, seed); + }}, +{ "Custom", [](unsigned int /*ignored_d*/, const BenchmarkConfig& cfg) { + if (cfg.custom_A_file.empty() || cfg.custom_b_file.empty()) { + throw std::runtime_error("\n[CRITICAL ERROR] Custom polytope chosen but CSV file paths are missing in the config JSON!\n"); + } + + // verify files exist on disk before touching them + if (!std::filesystem::exists(cfg.custom_A_file)) { + throw std::runtime_error( + "\n[CRITICAL ERROR] Cannot find A matrix file: " + cfg.custom_A_file + + "\n-> Hint: Ensure the CSV is in the same directory you are running the executable from, or use an absolute path in your config.\n" + ); + } + if (!std::filesystem::exists(cfg.custom_b_file)) { + throw std::runtime_error( + "\n[CRITICAL ERROR] Cannot find b vector file: " + cfg.custom_b_file + + "\n-> Hint: Ensure the CSV is in the same directory you are running the executable from, or use an absolute path in your config.\n" + ); + } + + return load_custom_polytope(cfg.custom_A_file, cfg.custom_b_file); + }} + }; + + auto it = factory.find(choice); + if (it != factory.end()) { + return it->second(dim, config); // Call the matched lambda + } else { + std::cerr << "!!! Unknown polytope choice: " << choice << ". Defaulting to Cube.\n"; + return generate_cube(dim, false); + } +} \ No newline at end of file diff --git a/benchmark/include/progress_bar.hpp b/benchmark/include/progress_bar.hpp new file mode 100644 index 000000000..b020790cb --- /dev/null +++ b/benchmark/include/progress_bar.hpp @@ -0,0 +1,99 @@ +#pragma once + +#include +#include +#include +#include +#include + +#if defined(_WIN32) +#include +#else +#include +#include +#endif + +inline unsigned int get_terminal_width(unsigned int fallback = 80) { +#if defined(_WIN32) + CONSOLE_SCREEN_BUFFER_INFO csbi; + if (GetConsoleScreenBufferInfo(GetStdHandle(STD_OUTPUT_HANDLE), &csbi)) { + return static_cast(csbi.srWindow.Right - csbi.srWindow.Left + 1); + } + return fallback; +#else + struct winsize w; + if (ioctl(STDOUT_FILENO, TIOCGWINSZ, &w) == 0 && w.ws_col > 0) { + return w.ws_col; + } + return fallback; +#endif +} + +inline void draw_progress_bar(const std::string& walk_name, + unsigned int current, + unsigned int total, + unsigned int total_samples_so_far, + unsigned int current_ESS, + unsigned int walk_len, + double elapsed_seconds, + double estimated_remaining_seconds, + int bar_width = 25) +{ + if (total == 0) return; + + float progress = static_cast(current) / total; + int pos = static_cast(bar_width * progress); + + std::ostringstream oss; + + oss << "[" << walk_name << "] Batch: ["; + for (int i = 0; i < bar_width; ++i) { + if (i < pos) oss << "="; + else if (i == pos) oss << ">"; + else oss << " "; + } + oss << "] " << static_cast(progress * 100.0) << "% (" + << current << "/" << total << ")"; + + if (current_ESS > 0) { + double live_mixing_ratio = static_cast(total_samples_so_far * walk_len) / current_ESS; + oss << " | Mix Ratio: " << std::fixed << std::setprecision(2) << live_mixing_ratio; + } else { + oss << " | Mix Ratio: N/A"; + } + + unsigned int e_total_secs = static_cast(elapsed_seconds); + unsigned int e_hours = e_total_secs / 3600; + unsigned int e_minutes = (e_total_secs % 3600) / 60; + unsigned int e_seconds = e_total_secs % 60; + + oss << " | Elapsed: "; + if (e_hours > 0) oss << e_hours << "h " << e_minutes << "m"; + else if (e_minutes > 0) oss << e_minutes << "m " << e_seconds << "s"; + else oss << e_seconds << "s"; + + if (estimated_remaining_seconds >= 0.0) { + unsigned int total_secs = static_cast(estimated_remaining_seconds); + unsigned int hours = total_secs / 3600; + unsigned int minutes = (total_secs % 3600) / 60; + unsigned int seconds = total_secs % 60; + + oss << " | ETA: "; + if (hours > 0) oss << hours << "h " << minutes << "m"; + else if (minutes > 0) oss << minutes << "m " << seconds << "s"; + else oss << seconds << "s"; + } else { + oss << " | ETA: Calculating..."; + } + + std::string line = oss.str(); + + unsigned int term_width = get_terminal_width(); + unsigned int max_width = (term_width > 1) ? term_width - 1 : term_width; + + if (line.size() > max_width) { + line = line.substr(0, max_width); + } + + std::cout << "\r" << std::string(term_width, ' ') << "\r" << line << std::flush; +} \ No newline at end of file diff --git a/benchmark/include/rounding.hpp b/benchmark/include/rounding.hpp new file mode 100644 index 000000000..878df3673 --- /dev/null +++ b/benchmark/include/rounding.hpp @@ -0,0 +1,61 @@ +#ifndef ROUNDING_HPP +#define ROUNDING_HPP + +#include +#include +#include +#include + +#include "core_types.hpp" +#include "inscribed_ellipsoid_rounding.hpp" + +template +void apply_polytope_rounding(const std::string& method, + PolytopeType& Polytope, + PointType& center, + MT& T, + VT& shift, + NT& round_val, + bool show_console_logs) { + + if (show_console_logs) { + std::cout << "[ROUNDING] Rounding is enabled. Applying " << method << " rounding...\n"; + } + + if (method == "max_ellipsoid") { + auto rounding_result = inscribed_ellipsoid_rounding(Polytope, center); + T = std::get<0>(rounding_result); + shift = std::get<1>(rounding_result); + round_val = std::get<2>(rounding_result); + } + else if (method == "log_barrier") { + auto rounding_result = inscribed_ellipsoid_rounding(Polytope, center); + T = std::get<0>(rounding_result); + shift = std::get<1>(rounding_result); + round_val = std::get<2>(rounding_result); + } + else if (method == "vaidya_barrier") { + auto rounding_result = inscribed_ellipsoid_rounding(Polytope, center); + T = std::get<0>(rounding_result); + shift = std::get<1>(rounding_result); + round_val = std::get<2>(rounding_result); + } + else if (method == "volumetric_barrier") { + auto rounding_result = inscribed_ellipsoid_rounding(Polytope, center); + T = std::get<0>(rounding_result); + shift = std::get<1>(rounding_result); + round_val = std::get<2>(rounding_result); + } + else { + throw std::runtime_error("Unknown rounding method: " + method); + } + + // Since rounding shifts the polytope to the origin we will use a zero vector as the center. + center = PointType(VT::Zero(Polytope.dimension())); + + if (show_console_logs) { + std::cout << "[ROUNDING] Rounding complete. Round value: " << round_val << "\n\n"; + } +} + +#endif // ROUNDING_HELPER_HPP \ No newline at end of file diff --git a/benchmark/include/walk_adapters.hpp b/benchmark/include/walk_adapters.hpp new file mode 100644 index 000000000..ee717991c --- /dev/null +++ b/benchmark/include/walk_adapters.hpp @@ -0,0 +1,342 @@ +#pragma once + +#include "core_types.hpp" +#include "benchmark_utils.hpp" +#include "walk_parameters.hpp" +#include "random_walks/random_walks.hpp" + +#include "sampling/sampling.hpp" +/* + * ----------------- + * This file provides a uniform interface for all random walk implementations + * used. + * + * Different walk algorithms expose different constructor signatures + * and apply() methods. For example, most uniform samplers can be constructed + * as Walk(P, p, rng), while Gaussian samplers require additional + * parameters (e.g. parameter a_i) and different apply() signatures. + * + * To avoid having special cases throughout the benchmark code, this file + * introduces the WalkAdapter abstraction. Each adapter exposes a common + * apply_batch() function that: + * + * 1. Constructs the corresponding walk object. + * 2. Executes a specified number of walk steps. + * 3. Collects generated sample points into a batch. + * 4. Enforces benchmark time limits. + * + * The overall philosophy is that every random walk method can be executed + * through the same adapter interface, allowing the rest of the framework + * (benchmark runners, statistics collection, batching logic, etc.) to remain + * completely independent of the underlying sampling algorithm. + * Check walk_run.hpp for the universal run function. + * MACROS are used in an attempt to shorten the size of this file since most of the code is the same anyway +*/ + +// Uniform Walk Type definitions +typedef BallWalk::template Walk BallWalkType; +typedef BilliardWalk::template Walk BilliardWalkType; +typedef AcceleratedBilliardWalk::template Walk AcceleratedBilliardWalkType; +typedef SparseBilliardWalk::template Walk SparseBilliardWalkType; +typedef CDHRWalk::template Walk CDHRWalkType; +typedef RDHRWalk::template Walk RDHRWalkType; +typedef DikinWalk::template Walk DikinWalkType; +typedef JohnWalk::template Walk JohnWalkType; +typedef VaidyaWalk::template Walk VaidyaWalkType; + +typedef GaussianBallWalk::template Walk GaussianBallWalkType; +typedef GaussianCDHRWalk::template Walk GaussianCDHRWalkType; + +typedef ShakeAndBakeWalk::template Walk ShakeAndBakeWalkType; +typedef BilliardShakeAndBakeWalk::template Walk BilliardSBWalkType; + +typedef BCDHRWalk::template Walk BCDHRWalkType; +typedef BRDHRWalk::template Walk BRDHRWalkType; + + +// For Uniform Walks +template +struct WalkAdapter { + + static constexpr bool supports_chunking = true; + + // Initialization + static WalkType init(HPOLYTOPE& P, Point& p, const BenchmarkConfig& config, RNGType& rng) { + return WalkType(P, p, rng); + } + + // apply batch + static void apply_batch(WalkType& walk, HPOLYTOPE& P, Point& p, unsigned int batch_size, + unsigned int walk_len, std::vector& batchPoints, + const BenchmarkConfig& config, RNGType& rng, Timer& walk_timer) + { + for (unsigned int i = 0; i < batch_size; ++i) { + walk.apply(P, p, walk_len, rng); + batchPoints.push_back(p); + + if (i % 50 == 0 && walk_timer.get_total_time() > config.time_limit_sec) { + break; + } + } + } +}; + +// For GAUSSIAN WALKs +#define REGISTER_GAUSSIAN_ADAPTER(WALK_TYPE, JSON_NAME) \ +template <> \ +struct WalkAdapter { \ + \ + static constexpr bool supports_chunking = true; \ + \ + /* Initialization */ \ + static WALK_TYPE init(HPOLYTOPE& P, Point& p, const BenchmarkConfig& config, RNGType& rng) { \ + double a_i = 1.0; \ + auto walk_iter = config.walk_settings.find(JSON_NAME); \ + if (walk_iter != config.walk_settings.end()) { \ + a_i = walk_iter->second.a_i_param; \ + } \ + return WALK_TYPE(P, p, a_i, rng); \ + } \ + \ + /* aplly batch */ \ + static void apply_batch(WALK_TYPE& walk, HPOLYTOPE& P, Point& p, unsigned int batch_size, \ + unsigned int walk_len, std::vector& batchPoints, \ + const BenchmarkConfig& config, RNGType& rng, \ + Timer& walk_timer) \ + { \ + /* We fetch a_i to pass into the apply method */ \ + double a_i = 1.0; \ + auto walk_iter = config.walk_settings.find(JSON_NAME); \ + if (walk_iter != config.walk_settings.end()) { \ + a_i = walk_iter->second.a_i_param; \ + } \ + \ + for (unsigned int i = 0; i < batch_size; ++i) { \ + walk.apply(P, p, a_i, walk_len, rng); \ + batchPoints.push_back(p); \ + \ + /* The Timeout Check */ \ + if (i % 50 == 0 && walk_timer.get_total_time() > config.time_limit_sec) { \ + break; \ + } \ + } \ + } \ +}; + +REGISTER_GAUSSIAN_ADAPTER(GaussianBallWalkType, "GaussianBallWalk") +REGISTER_GAUSSIAN_ADAPTER(GaussianCDHRWalkType, "GaussianCDHRWalk") + +// FOR Billiard Shake-and-Bake walks +// For these walks we need a point on a facet to start with. +// To that end, we use the billiard logic to shoot a ray from our interior initial point and check where we hit a facet +template <> +struct WalkAdapter { + + static constexpr bool supports_chunking = true; + + // Initialization + static BilliardSBWalkType init(HPOLYTOPE& P, Point& p, const BenchmarkConfig& config, RNGType& rng) { + int nr = 10; // Hardcoded number of reflections + + unsigned int n = P.dimension(); + Point v = GetDirection::apply(n, rng); + + // Temporary structures required by line_positive_intersect + typename Point::Coeff lambdas(P.num_of_hyperplanes()); + typename Point::Coeff Av(P.num_of_hyperplanes()); + lambdas.setZero(); + Av.setZero(); + + // Find intersection + std::pair pbpair = P.line_positive_intersect(p, v, lambdas, Av); + + // Permanently move the starting point to the boundary + p += (pbpair.first * v); + int initial_facet = pbpair.second; + + // Construct and return the walk + return BilliardSBWalkType(P, p, rng, initial_facet, nr); + } + + // apply batch + static void apply_batch( + BilliardSBWalkType& walk, + HPOLYTOPE& P, + Point& p, + unsigned int batch_size, + unsigned int walk_len, + std::vector& batchPoints, + const BenchmarkConfig& config, + RNGType& rng, + Timer& walk_timer) + { + for (unsigned int i = 0; i < batch_size; ++i) { + + // Advance the already-initialized walk + walk.apply(P, walk_len, rng); + + p = walk.getCurrentPoint(); + batchPoints.push_back(p); + + // Timeout Check + if (i % 50 == 0 && walk_timer.get_total_time() > config.time_limit_sec) { + break; + } + } + } +}; + +// FOR Shake-and-Bake walks +template <> +struct WalkAdapter { + + static constexpr bool supports_chunking = true; + + // Initialization. Find the boundary point + static ShakeAndBakeWalkType init(HPOLYTOPE& P, Point& p, const BenchmarkConfig& config, RNGType& rng) { + int initial_facet = 0; + auto b = P.get_vec(); + bool on_boundary = false; + + // Check if the current point 'p' is already on a facet + for (int i = 0; i < P.num_of_hyperplanes(); ++i) { + if (std::abs(P.get_row(i).dot(p.getCoefficients()) - b(i)) < 1e-7) { + initial_facet = i; + on_boundary = true; + break; + } + } + + // If it's an interior point, project it to the boundary + if (!on_boundary) { + typename Point::Coeff v_vec = Point::Coeff::Zero(P.dimension()); + v_vec(0) = 1.0; + + double min_lambda = std::numeric_limits::max(); + + for (int i = 0; i < P.num_of_hyperplanes(); ++i) { + double v_dot_a = P.get_row(i).dot(v_vec); + if (v_dot_a > 1e-10) { + double dist = b(i) - P.get_row(i).dot(p.getCoefficients()); + double lam = dist / v_dot_a; + + if (lam > 0 && lam < min_lambda) { + min_lambda = lam; + initial_facet = i; + } + } + } + + // Move the initial point permanently to the collision spot + typename Point::Coeff new_coords = p.getCoefficients() + (min_lambda * v_vec); + p = Point(new_coords); + } + + // Construct and return the walk + return ShakeAndBakeWalkType(P, p, initial_facet, rng); + } + + // Update the walk + static void apply_batch(ShakeAndBakeWalkType& walk, HPOLYTOPE& P, Point& p, + unsigned int batch_size, unsigned int walk_len, + std::vector& batchPoints, const BenchmarkConfig& config, + RNGType& rng, Timer& walk_timer) + { + for (unsigned int i = 0; i < batch_size; ++i) { + + walk.apply(P, walk_len, rng); + + p = walk.getCurrentPoint(); + batchPoints.push_back(p); + + if (i % 50 == 0 && walk_timer.get_total_time() > config.time_limit_sec) { + break; + } + } + } +}; + + +// FOR BOUNDARY HIT-AND-RUN WALKS (BCDHR, BRDHR) +#define REGISTER_BOUNDARY_HR_ADAPTER(WALK_TYPE, JSON_NAME) \ +template <> \ +struct WalkAdapter { \ + \ + static constexpr bool supports_chunking = true; \ + \ + /* Initialization */ \ + static WALK_TYPE init(HPOLYTOPE& P, Point& p, const BenchmarkConfig& config, RNGType& rng) { \ + /* Constructor with 3 arguments */ \ + return WALK_TYPE(P, p, rng); \ + } \ + \ + /* Apply batch */ \ + static void apply_batch(WALK_TYPE& walk, HPOLYTOPE& P, Point& p, unsigned int batch_size, \ + unsigned int walk_len, std::vector& batchPoints, \ + const BenchmarkConfig& config, RNGType& rng, \ + Timer& walk_timer) \ + { \ + /* Dummy points to catch the boundary chord endpoints for this batch */ \ + Point chord_p1 = p; \ + Point chord_p2 = p; \ + \ + for (unsigned int i = 0; i < batch_size; ++i) { \ + /* Advance the preserved walk object */ \ + walk.apply(P, chord_p1, chord_p2, walk_len, rng); \ + \ + /* Extract the updated internal point using getter */ \ + p = walk.getCurrentPoint(); \ + batchPoints.push_back(p); \ + \ + /* Timeout Check */ \ + if (i % 50 == 0 && walk_timer.get_total_time() > config.time_limit_sec) { \ + break; \ + } \ + } \ + } \ +}; + +REGISTER_BOUNDARY_HR_ADAPTER(BCDHRWalkType, "BCDHRWalk") +REGISTER_BOUNDARY_HR_ADAPTER(BRDHRWalkType, "BRDHRWalk") + +// Riemannian Hamiltonian +// This method works a bit differently so we run a big chunk of points at once. +template <> +struct WalkAdapter { + + static constexpr bool supports_chunking = false; + + // Dummy initialization so that we satisfy the generic template caller + static int init(HPOLYTOPE& P, Point& p, const BenchmarkConfig& config, RNGType& rng) { + return 0; + } + + // apply_batch takes the dummy integer, which we safely ignore + static void apply_batch(int& dummy_state, HPOLYTOPE& P, Point& p, unsigned int batch_size, + unsigned int walk_len, std::vector& batchPoints, + const BenchmarkConfig& config, RNGType& rng, Timer& walk_timer) + { + using Func = ZeroScalarFunctor; + using Grad = ZeroFunctor; + using Hess = ZeroFunctor; + + Func f; + Grad g; + Hess h; + + std::list temp_list; + int n_burns = 1000; + + // execute_crhmc handles the massive chunk and burn-in all at once + execute_crhmc, Grad, Func, Hess, CRHMCWalk, 1>( + P, rng, temp_list, 1, batch_size, n_burns, &g, &f, &h + ); + + batchPoints.insert(batchPoints.end(), temp_list.begin(), temp_list.end()); + + if (!temp_list.empty()) { + p = temp_list.back(); + } + } +}; + diff --git a/benchmark/include/walk_parameters.hpp b/benchmark/include/walk_parameters.hpp new file mode 100644 index 000000000..aefb71dfe --- /dev/null +++ b/benchmark/include/walk_parameters.hpp @@ -0,0 +1,51 @@ +#pragma once + +#include +#include +#include + +// In this file you can find functions for configuration structures and helper +// functions used to control benchmark experiments for random walks. + +// A struct to hold the specific settings for a single random walk +struct WalkSettings { + bool enabled; + unsigned int samples; + unsigned int walk_len_multiplier; + unsigned int walk_len_base; + + // For Gaussian + double a_i_param = 1.0; +}; + +// The global configuration object +struct BenchmarkConfig { + + unsigned int target_ESS; + double time_limit_sec; + int base_seed; + unsigned int dimension; + std::vector dimensions; + double angle; + std::string polytope_choice; + std::string custom_A_file; + std::string custom_b_file; + bool use_dynamic_batch; + bool write_to_file; + bool rounding; + std::string rounding_method; + bool auto_walk; + bool show_console_logs; + bool show_menu; + + std::map walk_settings; +}; + + + +// Function declaration to load the JSON file +BenchmarkConfig load_benchmark_config(const std::string& filepath); + +// Helper functions that the runners will use +unsigned int get_initial_batch_size(const std::string& walk_name, const BenchmarkConfig& config); +unsigned int compute_dynamic_walk_len(const std::string& walk_name, unsigned int dim, const BenchmarkConfig& config); \ No newline at end of file diff --git a/benchmark/include/walk_registry.hpp b/benchmark/include/walk_registry.hpp new file mode 100644 index 000000000..628159a19 --- /dev/null +++ b/benchmark/include/walk_registry.hpp @@ -0,0 +1,69 @@ +#pragma once + +#include +#include +#include + +#include "core_types.hpp" +#include "walk_run.hpp" +#include "walk_result.hpp" + +/* + * This file allows sampling methods to be selected by name at runtime, + * instead of creating specific walk types directly in the code. + * + * Every walk is registered with a string name and exposed through the same + * function interface. This makes it possible to store and execute different + * walk implementations in a uniform way. + * + * RunFunction defines the common function signature that all registered walks + * must follow. + * + * WalkRegistry stores the mapping between walk names and their corresponding + * execution functions. It acts as a lookup table, making it easy to choose a + * sampling method based on user input or benchmark settings. + * + * The execute_walk() helper adapts template-based walk implementations so they + * can be stored in the registry and called through the common interface. + * + * Main functions: + * + * - register_walk(): Registers a walk in the registry. + * - get_walk_registry(): Returns the global registry. + * - initialize_all_walks(): Registers all available walks. + * + * Together with walk_adapters.hpp, this file provides the infrastructure for + * selecting and running sampling methods through a single, consistent interface. + */ + +// Each walk must match this callable signature +using RunFunction = std::function; + +// Registry type alias +using WalkRegistry = std::map; + +// Function that returns the global registry +WalkRegistry& get_walk_registry(); + +// Helper to register a walk. It binds the name to the specific function. +void register_walk(const std::string& name, RunFunction fn); + +// Executes the walk by calling the sample_using_walk +template +WalkResult execute_walk( + HPOLYTOPE& P, + const Point& c, + RNGType& r, + const BenchmarkConfig& cfg, + const std::string& name +) { + return sample_using_walk(P, c, r, cfg, name); +} + +void initialize_all_walks(); \ No newline at end of file diff --git a/benchmark/include/walk_result.hpp b/benchmark/include/walk_result.hpp new file mode 100644 index 000000000..6a1ee0c83 --- /dev/null +++ b/benchmark/include/walk_result.hpp @@ -0,0 +1,39 @@ +#pragma once + +#include "core_types.hpp" +#include +#include + +// In this file you can find the structures and functions responsible for +// post-processing and summarizing the results of random walk sampling. + +struct WalkResult { + std::vector samples; + unsigned int final_ess; + double generation_time; + double ess_time; + unsigned int walk_len; +}; + + +struct WalkStatistics { + unsigned int final_ess; + double max_psrf; + double ks_statistic; + double ks_p_value; + double total_time; + double ess_time; + double mixing_ratio; +}; + +WalkStatistics process_and_print_results( + const std::vector& samples, + HPOLYTOPE& polytope, + const std::string& walk_name, + double total_generation_time, + unsigned int precalculated_ess, + double total_ess_time, + unsigned int walk_len, + const std::string& polytope_name, + bool show_console_logs +); \ No newline at end of file diff --git a/benchmark/include/walk_run.hpp b/benchmark/include/walk_run.hpp new file mode 100644 index 000000000..4f9b3684d --- /dev/null +++ b/benchmark/include/walk_run.hpp @@ -0,0 +1,216 @@ +#pragma once + +#include "core_types.hpp" +#include "walk_parameters.hpp" +#include "benchmark_utils.hpp" +#include "walk_result.hpp" +#include "walk_adapters.hpp" +#include "diagnostics.hpp" +#include "progress_bar.hpp" +#include "dynamic_batch_size.hpp" + +#include "sampling/random_point_generators.hpp" +#include +#include + +// In this file you can find sample using walk function used to generate +// points from a polytope using a specified random walk algorithm. + +template +WalkResult sample_using_walk(HPOLYTOPE& Polytope, + Point const& start_point, + RNGType& rng, + const BenchmarkConfig& config, + const std::string& walk_name) +{ + // Initial point and dimension + Point starting_point = start_point; + unsigned int dim = Polytope.dimension(); + + // ESS placeholder and loop number + unsigned int current_ESS = 0; + unsigned int loop_step = 1; + unsigned int previous_ESS = 0; + double estimated_remaining_seconds = -1.0; + + // Random generator and placeholder matrix for samples + typedef RandomPointGenerator Generator; + std::vector allSamples; + allSamples.reserve(config.target_ESS * 10); + + // Find initial batch size and walk_len based on config file and user choices + unsigned int batch_size = get_initial_batch_size(walk_name, config); + unsigned int walk_len = compute_dynamic_walk_len(walk_name, dim, config); + + // initiallize the walk + auto walk = WalkAdapter::init(Polytope, starting_point, config, rng); + + // Timer objects that stored the time + Timer walk_timer(walk_name); + Timer ess_timer(walk_name); + walk_timer.start(); + Timer global_timer("TotalTime"); + global_timer.start(); + + // Main while loop. We sample intil we hit target ESS. + while (current_ESS < config.target_ESS && global_timer.get_total_time() < config.time_limit_sec) { + + std::vector batchPoints; + + // here we calculate the live time estimate for the upcoming batch + double elapsed_walk_time = walk_timer.get_total_time(); + unsigned int samples_so_far = allSamples.size(); + + if (samples_so_far > 0 && elapsed_walk_time > 0.01 && current_ESS > 0) { + // how many samples we generate per second + double samples_per_sec = static_cast(samples_so_far) / elapsed_walk_time; + // how many samples are still needed to reach Target ESS + double remaining_ess_needed = static_cast(config.target_ESS) - current_ESS; + // Samples needed = (Remaining ESS) * (Average samples required per 1 ESS) + double estimated_remaining_samples = remaining_ess_needed * (static_cast(samples_so_far) / current_ESS); + // ETA = Remaining Samples / Sample Rate + estimated_remaining_seconds = estimated_remaining_samples / samples_per_sec; + } else { + estimated_remaining_seconds = -1.0; // Still on the first batch / no ESS data yet + } + + // Progress Bar + unsigned int chunk_size = 250; // How many samples to generate before updating the bar + unsigned int generated_this_batch = 0; + + if constexpr (WalkAdapter::supports_chunking) { + + unsigned int chunk_size = 250; + unsigned int generated_this_batch = 0; + unsigned int samples_before_this_batch = allSamples.size(); + + while (generated_this_batch < batch_size) { + unsigned int current_chunk = std::min(chunk_size, batch_size - generated_this_batch); + std::vector chunkPoints; + + // The main call. The same call is valid for all methods + WalkAdapter::apply_batch( + walk, Polytope, starting_point, current_chunk, walk_len, chunkPoints, config, rng, walk_timer + ); + + if (!chunkPoints.empty()) { + starting_point = chunkPoints.back(); + } + allSamples.insert(allSamples.end(), chunkPoints.begin(), chunkPoints.end()); + generated_this_batch += chunkPoints.size(); + double global_elapsed_seconds = global_timer.get_total_time(); + + // Draw progress bar with mixing ratio + draw_progress_bar( + walk_name, + generated_this_batch, + batch_size, + samples_before_this_batch, + current_ESS, + walk_len, + global_elapsed_seconds, + estimated_remaining_seconds + ); + + if (walk_timer.get_total_time() > config.time_limit_sec) { + break; + } + } + } else { + // RIEMANNIAN METHOD: Massive Single Batch + std::cout << "\r" << std::string(100, ' ') << "\r[" << walk_name + << "] Generating massive batch of " << batch_size + << " points (Tuning physics engine)..." << std::flush; + + if (estimated_remaining_seconds >= 0.0) { + std::cout << " (ETA: " << static_cast(estimated_remaining_seconds) << "s remaining)..." << std::flush; + } + + std::vector singleBatchPoints; + WalkAdapter::apply_batch( + walk, Polytope, starting_point, batch_size, walk_len, singleBatchPoints, config, rng, walk_timer + ); + + if (!singleBatchPoints.empty()) { + starting_point = singleBatchPoints.back(); + } + allSamples.insert(allSamples.end(), singleBatchPoints.begin(), singleBatchPoints.end()); + + if (walk_timer.get_total_time() > config.time_limit_sec) { + std::cout << "\n[" << walk_name << "] Time limit exceeded during batch generation.\n"; + break; + } + } + + std::cout << "\r" << std::string(100, ' ') << "\r[" << walk_name << "] Batch Complete! Calculating ESS..." << std::flush; + + walk_timer.stop(""); + ess_timer.start(); + + // Calculate ESS + MT samples = vector_to_eigen(allSamples); + current_ESS = compute_ess(samples); + + ess_timer.stop(""); + walk_timer.start(); + + std::cout << "\r" << std::string(120, ' ') << "\r[" << walk_name << "] Samples: " << allSamples.size() + << " | ESS: " << current_ESS; + + // If the user has dynamic_batch_size off we end here after "samples" number of samples are generated. + if (!config.use_dynamic_batch) { + std::cout << "\n[" << walk_name << "] Generated " << allSamples.size() + << " fixed samples. Stopping without checking ESS.\n"; + break; + } + + // Main check. If we passed the target stop immediatelly. + if (current_ESS >= config.target_ESS) { + std::cout << "\n[" << walk_name << "] Reached Target ESS. Stopping.\n"; + break; + } + + if (global_timer.get_total_time() > config.time_limit_sec) { + std::cout << "\n[" << walk_name << "] Time limit exceeded. Stopping.\n"; + break; + } + + // THE DYNAMIC BATCH SIZE + if (config.use_dynamic_batch) { + double remaining_time_sec = config.time_limit_sec - global_timer.get_total_time(); + double current_samples_per_sec = (elapsed_walk_time > 0.01) + ? static_cast(allSamples.size()) / elapsed_walk_time + : 0.0; + + batch_size = compute_next_batch_size( + config.target_ESS, + current_ESS, + previous_ESS, + allSamples.size(), + batch_size, + config.dimension, + remaining_time_sec, + current_samples_per_sec + ); + } + + std::cout << " | Next Batch: " << batch_size << std::flush; + + previous_ESS = current_ESS; + + // Failsafe + if (loop_step > 30 && current_ESS < config.target_ESS) { + std::cout << "\n[" << walk_name << "] I am sorry. I did not converge!\n"; + break; + } + + loop_step++; + } + + std::cout << "\n[" << walk_name << "] DONE. Samples Generated: " << allSamples.size() << "\n"; + //global_timer.stop("Total Execution Time"); + double final_gen_time = walk_timer.get_total_time(); + double final_ess_time = ess_timer.get_total_time(); + + return { allSamples, current_ESS, final_gen_time, final_ess_time, walk_len }; +} diff --git a/benchmark/src/benchmark_main.cpp b/benchmark/src/benchmark_main.cpp new file mode 100644 index 000000000..e34ea77f8 --- /dev/null +++ b/benchmark/src/benchmark_main.cpp @@ -0,0 +1,5 @@ +#include "benchmark_run.hpp" + +int main(int argc, char** argv) { + return run_benchmark(argc, argv); +} \ No newline at end of file diff --git a/benchmark/src/benchmark_run.cpp b/benchmark/src/benchmark_run.cpp new file mode 100644 index 000000000..6ea8b914b --- /dev/null +++ b/benchmark/src/benchmark_run.cpp @@ -0,0 +1,257 @@ +#include +#include +#include +#include + +#include "../include/core_types.hpp" +#include "../include/walk_parameters.hpp" +#include "../include/benchmark_utils.hpp" +#include "../include/geometry_utils.hpp" +#include "../include/walk_result.hpp" +#include "../include/walk_registry.hpp" +#include "../include/menu.hpp" +#include "../include/rounding.hpp" + +#include "../include/polytope_generation.hpp" +#include "known_polytope_generators.h" +#include "order_polytope_generator.h" + +#include "benchmark_run.hpp" + +using namespace std; + +namespace po = boost::program_options; + +int run_benchmark(int argc, char** argv) { + + initialize_all_walks(); + + string config_file = "../config/walk_config.json"; // Base default + unsigned int dimension; + string walk_choice; + string polytope_cli_choice; + + // We first pre-parse only the config file flag so we can load the JSON defaults first + for (int i = 1; i < argc; ++i) { + string arg = argv[i]; + if ((arg == "-c" || arg == "--config") && i + 1 < argc) { + config_file = argv[i + 1]; + break; + } + } + + // Load the Configuration + cout << "Loading configuration from: " << config_file << "\n"; + BenchmarkConfig config = load_benchmark_config(config_file); + + // Setup Command Line Options (using JSON values as our defaults) + po::options_description desc("Benchmark Options"); + desc.add_options() + ("help,h", "Produce help message") + ("config,c", po::value(&config_file)->default_value(config_file), "Path to JSON config") + ("dim,d", po::value(&dimension)->default_value(config.dimension), "Dimension of the polytope") + ("polytope,p", po::value(&polytope_cli_choice)->default_value(config.polytope_choice), "Choose polytope: Cube, Simplex, Birkhoff, Cross, OrderPolytope, Custom") + ("walk,w", po::value(&walk_choice)->default_value("All"), + "Specific walk to run, or 'All'. Valid options:\n" + " - BallWalk\n" + " - BilliardWalk\n" + " - AcceleratedBilliardWalk\n" + " - SparseBilliardWalk\n" + " - CDHRWalk\n" + " - RDHRWalk\n" + " - DikinWalk\n" + " - JohnWalk\n" + " - VaidyaWalk\n" + " - GaussianBallWalk\n" + " - GaussianCDHRWalk\n" + " - BilliardShakeAndBakeWalk\n" + " - ShakeAndBakeWalk\n" + " - BCDHRWalk\n" + " - BRDHRWalk\n" + " - CRHMCWalk"); + + // Error-Checking Block + po::variables_map vm; + try { + po::store(po::parse_command_line(argc, argv, desc), vm); + po::notify(vm); + } catch (const exception& e) { + cerr << "Error parsing arguments: " << e.what() << "\n"; + return 1; + } + + if (vm.count("help")) { + cout << desc << "\n"; + return 0; + } + + config.polytope_choice = polytope_cli_choice; + + // --- NEW MENU LOGIC --- + if (config.show_menu) { + bool continue_to_benchmark = run_interactive_menu(config, walk_choice); + if (!continue_to_benchmark) { + return 0; // Exit gracefully if they chose 3 + } + } + // --------------------------- + + std::vector dimensions_to_run; + if (config.polytope_choice == "Custom") { + dimensions_to_run = { 0 }; + } else { + dimensions_to_run = config.dimensions; + } + + for (unsigned int current_dim : dimensions_to_run) { + + config.dimension = current_dim; + HPOLYTOPE Polytope_simple; + try { + Polytope_simple = create_polytope(config.polytope_choice, config.dimension, config); + } + catch (const std::exception& e) { + std::cerr << e.what() << "\n"; + std::cerr << ">>> Error. Please fix the config or file paths.\n"; + return 1; + } + + config.dimension = Polytope_simple.dimension(); + dimension = config.dimension; + + double angle = config.angle; + + // Print basic info + if (config.show_console_logs) { + cout << "Target ESS: " << config.target_ESS << "\n"; + cout << "Dimension: " << dimension << "\n"; + cout << "Polytope: " << config.polytope_choice << "\n"; + cout << "Rotation angle is: " << angle << "\n"; + cout << "Dynamic batch size is on: " << config.use_dynamic_batch << "\n"; + cout << "Rounding is on: " << config.rounding << "\n"; + cout << "Auto-walk is on: " << (config.auto_walk ? "true" : "false") << "\n"; + cout << "\n" << string(40, '=') << "\n"; + cout << "*** Running for dimension " << dimension << " ***\n"; + } + // We need this copy to pass the original polytope to the metrics if roundeing was on. + HPOLYTOPE Polytope_rotated = rotate_all_dims(Polytope_simple, angle); + HPOLYTOPE Polytope = Polytope_rotated; + + auto inner = Polytope.ComputeInnerBall(); + Point center = inner.first; + + // Variables to store transformation data + MT T; + VT shift; + NT round_val = 1.0; + + // ****ROUNDING***** + if (config.rounding) { + apply_polytope_rounding( + config.rounding_method, + Polytope, + center, + T, + shift, + round_val, + config.show_console_logs + ); + } + // -------------------------- + + // Setup RNG + RNGType rng(Polytope.dimension()); + + auto run_method = [&](const string& method_name) { + + // Access registry + auto& registry = get_walk_registry(); + auto walk_it = registry.find(method_name); + + if (walk_it != registry.end()) { + + WalkResult result = walk_it->second( + Polytope, + center, + rng, + config, + method_name + ); + + if (!result.samples.empty()) { + + // Reverse rounding + if (config.rounding) { + for (auto& pt : result.samples) { + + // Apply the reverse transformation: T * vector + shift + pt = T * pt.getCoefficients() + shift; + } + } + // ---------------------------------- + + // Write samples to txt file for later use + if (config.write_to_file) { + std::filesystem::create_directory("results"); + std::string filename = "results/" + config.polytope_choice + "_" + + std::to_string(dimension) + "_" + + method_name + "_samples.txt"; + + if (config.show_console_logs) { + std::cout << "[" << method_name << "] Saving " << result.samples.size() + << " points to " << filename << "...\n"; + } + + write_to_file(filename, result.samples); + if (config.show_console_logs) { + std::cout << "[" << method_name << "] File saved successfully.\n"; + } + } + + // Process results + process_and_print_results( + result.samples, + Polytope_rotated, + method_name, + result.generation_time, + result.final_ess, + result.ess_time, + result.walk_len, + config.polytope_choice, + config.show_console_logs + ); + } else { + cout << "!!! " << method_name << " failed to generate points.\n"; + } + + } else { + cout << "!!! Unknown walk type skipped: " << method_name << "\n"; + return; + } + }; + // Auto-walk handles the selection if enabled + if (config.auto_walk) { + string auto_selected_walk = determine_auto_walk(dimension); + if (config.show_console_logs) { + cout << "\n[Auto-Walk] Dimension " << dimension << " overriding config to run: " << auto_selected_walk << "\n"; + } + run_method(auto_selected_walk); + } + // If the user picked "All", iterate through the JSON keys. Otherwise, just run the one they requested. + else if (walk_choice == "All") { + for (const auto& pair : config.walk_settings) { + + if(pair.second.enabled) { + run_method(pair.first); + } + else { + // cout << "--- Skipping " << pair.first << " (Disabled in JSON) ---\n"; + } + } + } else { + run_method(walk_choice); + } + } + cout << "\nBenchmark Complete.\n"; + return 0; +} \ No newline at end of file diff --git a/benchmark/src/benchmark_utils.cpp b/benchmark/src/benchmark_utils.cpp new file mode 100644 index 000000000..13322c02e --- /dev/null +++ b/benchmark/src/benchmark_utils.cpp @@ -0,0 +1,69 @@ +#include "../include/benchmark_utils.hpp" +#include +#include + +PushBackWalkPolicy push_back_policy; + +// Timer Class. This is usefull to avoid using new time variables every time we need to count the time. +Timer::Timer(const std::string& name) : walk_name(name), total_time(0.0) {} + +void Timer::start() { + start_time = std::chrono::steady_clock::now(); + is_running = true; +} + +double Timer::stop(const std::string& label) { + auto end_time = std::chrono::steady_clock::now(); + double elapsed = std::chrono::duration(end_time - start_time).count(); + total_time += elapsed; + is_running = false; + + // Only print if a label was provided + if (!label.empty()) { + std::cout << "[" << walk_name << "] " << label << " = " << elapsed << " s\n"; + } + return elapsed; +} + +double Timer::get_total_time() const { + + if (is_running) { + auto current_time = std::chrono::steady_clock::now(); + double current_elapsed = std::chrono::duration(current_time - start_time).count(); + return total_time + current_elapsed; + } + + return total_time; +} + +void write_to_file(std::string filename, std::vector const& randPoints) { + std::ofstream out(filename); + if (!out.is_open()) { + std::cerr << "Error: Could not open " << filename << " for writing.\n"; + return; + } + + // Save current cout buffer and redirect to the file + auto coutbuf = std::cout.rdbuf(out.rdbuf()); + + for(size_t i = 0; i < randPoints.size(); ++i) { + randPoints[i].print(); + } + + // Reset cout back to standard output + std::cout.rdbuf(coutbuf); +} + +std::string determine_auto_walk(unsigned int dim) { + if (dim >= 1 && dim <= 10) { + return "BallWalk"; + } else if (dim >= 11 && dim <= 20) { + return "RDHRWalk"; + } else if (dim >= 21 && dim <= 30) { + return "BilliardWalk"; + } else if (dim >= 31 && dim <= 49) { + return "CDHRWalk"; + } else { + return "AcceleratedBilliardWalk"; + } +} \ No newline at end of file diff --git a/benchmark/src/walk_parameters.cpp b/benchmark/src/walk_parameters.cpp new file mode 100644 index 000000000..3cc7cd38f --- /dev/null +++ b/benchmark/src/walk_parameters.cpp @@ -0,0 +1,79 @@ +#include "../include/walk_parameters.hpp" +#include +#include +#include + +BenchmarkConfig load_benchmark_config(const std::string& filepath) { + BenchmarkConfig config; + boost::property_tree::ptree pt; + + try { + // Parse the JSON file + boost::property_tree::read_json(filepath, pt); + + // Read Global Settings + config.target_ESS = pt.get("global_settings.target_ESS", 500); + config.time_limit_sec = pt.get("global_settings.time_limit_sec", 1200.0); + config.base_seed = pt.get("global_settings.base_seed", 42); + + //dimensions array + auto dims_node = pt.get_child_optional("global_settings.dimensions"); + if (dims_node) { + for (const auto& item : *dims_node) { + config.dimensions.push_back(item.second.get_value()); + } + } else { + // fallback for older configs + config.dimensions.push_back(pt.get("global_settings.dimensions", 100)); + } + + config.polytope_choice = pt.get("global_settings.polytope_choice", "Cube"); + config.custom_A_file = pt.get("global_settings.custom_A_file", ""); + config.custom_b_file = pt.get("global_settings.custom_b_file", ""); + config.angle = pt.get("global_settings.rotation_angle", 0); + config.use_dynamic_batch = pt.get("global_settings.dynamic_batch_size", true); + config.write_to_file = pt.get("global_settings.write_to_file", false); + config.rounding = pt.get("global_settings.rounding", false); + config.rounding_method = pt.get("global_settings.rounding_method", "max_ellipsoid"); + config.auto_walk = pt.get("global_settings.auto_walk", false); + config.show_console_logs = pt.get("global_settings.show_console_logs", true); + config.show_menu = pt.get("global_settings.show_menu", false); + + // Iterate through the walks JSON object + for (const auto& walk_node : pt.get_child("walks")) { + std::string walk_name = walk_node.first; + WalkSettings settings; + + settings.enabled = walk_node.second.get("enabled", true); + settings.samples = walk_node.second.get("samples", 1000); + settings.walk_len_multiplier = walk_node.second.get("walk_len_multiplier", 0); + settings.walk_len_base = walk_node.second.get("walk_len_base", 1); + + //For Gaussian walks + settings.a_i_param = walk_node.second.get("a_i_param", 1.0); + + config.walk_settings[walk_name] = settings; + } + } catch (const boost::property_tree::ptree_error& e) { + std::cerr << "Error reading config file: " << e.what() << "\n"; + std::cerr << "Falling back to hardcoded defaults.\n"; + } + + return config; +} + +unsigned int get_initial_batch_size(const std::string& walk_name, const BenchmarkConfig& config) { + auto it = config.walk_settings.find(walk_name); + if (it != config.walk_settings.end()) { + return it->second.samples + (2 * config.dimension); + } + return 1000; +} + +unsigned int compute_dynamic_walk_len(const std::string& walk_name, unsigned int dim, const BenchmarkConfig& config) { + auto it = config.walk_settings.find(walk_name); + if (it != config.walk_settings.end()) { + return (dim * it->second.walk_len_multiplier) + it->second.walk_len_base; + } + return 1; +} \ No newline at end of file diff --git a/benchmark/src/walk_registry.cpp b/benchmark/src/walk_registry.cpp new file mode 100644 index 000000000..8d5b0a8ea --- /dev/null +++ b/benchmark/src/walk_registry.cpp @@ -0,0 +1,42 @@ +#include "walk_registry.hpp" +#include "walk_adapters.hpp" + +using namespace std; + +// Existing registry functions +WalkRegistry& get_walk_registry() { + static WalkRegistry registry; + return registry; +} + +void register_walk(const string& name, RunFunction fn) { + get_walk_registry()[name] = move(fn); +} + +void initialize_all_walks() { + // Uniform Walks + register_walk("BallWalk", execute_walk); + register_walk("BilliardWalk", execute_walk); + register_walk("AcceleratedBilliardWalk", execute_walk); + register_walk("SparseBilliardWalk", execute_walk); + register_walk("CDHRWalk", execute_walk); + register_walk("RDHRWalk", execute_walk); + register_walk("DikinWalk", execute_walk); + register_walk("JohnWalk", execute_walk); + register_walk("VaidyaWalk", execute_walk); + + // Gaussian Walks + register_walk("GaussianBallWalk", execute_walk); + register_walk("GaussianCDHRWalk", execute_walk); + + // Shake-and-Bake Walks + register_walk("ShakeAndBakeWalk", execute_walk); + register_walk("BilliardShakeAndBakeWalk", execute_walk); + + // Boundary Walks + register_walk("BCDHRWalk", execute_walk); + register_walk("BRDHRWalk", execute_walk); + + //Riemannian Walk + register_walk("CRHMCWalk", execute_walk); +} \ No newline at end of file diff --git a/benchmark/src/walk_result.cpp b/benchmark/src/walk_result.cpp new file mode 100644 index 000000000..6d92229f0 --- /dev/null +++ b/benchmark/src/walk_result.cpp @@ -0,0 +1,113 @@ +#include "../include/walk_result.hpp" +#include "../include/diagnostics.hpp" + +#include +#include +#include +#include +#include + +WalkStatistics process_and_print_results( + const std::vector& samples, + HPOLYTOPE& polytope, + const std::string& walk_name, + double total_generation_time, + unsigned int precalculated_ess, + double total_ess_time, + unsigned int walk_len, + const std::string& polytope_name, + bool show_console_logs) +{ + // Calculate mixing ratio (Total Steps / ESS) + double mixing_ratio = 0.0; + if (precalculated_ess > 0) { + // total steps = number of saved samples * thinning factor (walk_len) (no burn in included) + mixing_ratio = static_cast(samples.size() * walk_len) / precalculated_ess; + } + + // Calculate metrics + WalkStatistics stats = {precalculated_ess, 0.0, -1.0, -1.0, total_generation_time, total_ess_time, mixing_ratio}; + + if (samples.empty()) { + std::cerr << "[" << walk_name << "] Error: No samples to process!\n"; + return stats; + } + + //std::cout << "\n--- Processing Statistics for " << walk_name << " ---\n"; + + // Convert to Eigen Matrix + MT samples_mat = vector_to_eigen(samples); + + std::cout << std::fixed << std::setprecision(4); + + if (show_console_logs) { + // ESS and mixing rate + std::cout << "[" << walk_name << "] Final ESS: " << stats.final_ess << "\n"; + std::cout << "[" << walk_name << "] Mixing Ratio (Steps/ESS): " << stats.mixing_ratio << "\n"; + + // Time + std::cout << "[" << walk_name << "] Total Algorithm Time: " << stats.total_time << " seconds\n"; + std::cout << "[" << walk_name << "] Total ESS Time: " << stats.ess_time << " seconds\n"; + } + // PSRF + + stats.max_psrf = compute_psrf(samples); + + if (show_console_logs) { + std::cout << "[" << walk_name << "] Max PSRF: " << stats.max_psrf << "\n"; + } + + // Condition for KS Test + // Check if "Gaussian" is in the walk name because KS test is only for uniform + bool is_gaussian = (walk_name.find("Gaussian") != std::string::npos); + + if (!is_gaussian) { + auto ks_results = compute_ks_test(polytope, samples_mat, stats.final_ess); + stats.ks_statistic = ks_results.ks_stat; + stats.ks_p_value = ks_results.p_val; + + if (show_console_logs) { + std::cout << "[" << walk_name << "] KS Statistic: " << stats.ks_statistic << "\n"; + std::cout << "[" << walk_name << "] P-Value: " << stats.ks_p_value << "\n"; + } + } else { + if (show_console_logs) { + std::cout << "[" << walk_name << "] KS Test: Skipped (Gaussian Distribution)\n"; + } + } + std::cout << "--------------------------------------------------\n"; + + // Append results to the benchmark CSV file + std::string filename = "benchmark_results.csv"; + + // Check if file exists before we open it in append mode + bool file_exists = std::filesystem::exists(filename); + + std::ofstream outfile; + outfile.open(filename, std::ios_base::app); // Append mode + + if (outfile.is_open()) { + // If this is the very first time creating the file, write the header row + if (!file_exists) { + outfile << "Polytope,Dimension,Method,Time_Sec,Points,ESS,Mixing_Ratio,Max_PSRF,KS_Stat,KS_P_Value\n"; + } + + // Write the data row + outfile << polytope_name << ", " + << polytope.dimension() << ", " + << walk_name << ", " + << std::fixed << std::setprecision(4) << total_generation_time << ", " + << samples.size() << ", " + << stats.final_ess << ", " + << std::fixed << std::setprecision(4) << stats.mixing_ratio << ", " + << std::fixed << std::setprecision(4) << stats.max_psrf << ", " + << std::fixed << std::setprecision(4) << stats.ks_statistic << ", " + << std::fixed << std::setprecision(4) << stats.ks_p_value << "\n"; + + outfile.close(); + } else { + std::cerr << "!!! Unable to open " << filename << " to save data.\n"; + } + + return stats; +} \ No newline at end of file diff --git a/examples/general_sampling/.gitignore b/examples/general_sampling/.gitignore new file mode 100644 index 000000000..32f6d0bff --- /dev/null +++ b/examples/general_sampling/.gitignore @@ -0,0 +1,2 @@ +sampler +*walk.txt diff --git a/examples/general_sampling/CMakeLists.txt b/examples/general_sampling/CMakeLists.txt new file mode 100644 index 000000000..defcc8384 --- /dev/null +++ b/examples/general_sampling/CMakeLists.txt @@ -0,0 +1,66 @@ +cmake_minimum_required(VERSION 3.11) +project(VolEstiGeneralSampling) + +# -------------------------------------------------------- +# 1. Options & Flags +# -------------------------------------------------------- +option(DISABLE_NLP_ORACLES "Disable non-linear oracles" ON) +option(BUILTIN_EIGEN "Use eigen from ../../external" OFF) + +set(CMAKE_CXX_STANDARD 17) +set(CMAKE_CXX_STANDARD_REQUIRED ON) + +# Optimization & Standard Flags +add_compile_options(-O3 -march=native -DTIME_KEEPING -DBOOST_NO_AUTO_PTR -DEIGEN_NO_DEBUG) + +if(DISABLE_NLP_ORACLES) + add_definitions(-DDISABLE_NLP_ORACLES) +endif() + +# -------------------------------------------------------- +# 2. Dependencies (Using Relative Paths) +# -------------------------------------------------------- +# We use the helper scripts located in the main library +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() # This sets the variable ${QD_LIB} + +# -------------------------------------------------------- +# 3. Find Libraries Manually +# -------------------------------------------------------- +# Find lp_solve (standard paths + typical linux paths) +find_library(LP_SOLVE NAMES liblpsolve55.so lpsolve55 PATHS + /usr/lib/lp_solve + /usr/lib/x86_64-linux-gnu + /usr/local/lib +) + +if (NOT LP_SOLVE) + message(FATAL_ERROR "lp_solve library not found! Please install it.") +else () + message(STATUS "Library lp_solve found: ${LP_SOLVE}") +endif() + +# -------------------------------------------------------- +# 4. Include Directories +# -------------------------------------------------------- +# This gives you access to "volesti/include/..." +include_directories(BEFORE ../../include) +include_directories(BEFORE ../../external) +include_directories(BEFORE ../../external/minimum_ellipsoid) +include_directories(BEFORE ../../include/convex_bodies) +include_directories(BEFORE ../../include/generators) + +# -------------------------------------------------------- +# 5. Targets +# -------------------------------------------------------- +add_executable(general_sampling_time_limit general_sampling_time_limit.cpp) +target_link_libraries(general_sampling_time_limit ${LP_SOLVE} QD_LIB pthread) diff --git a/examples/general_sampling/README.md b/examples/general_sampling/README.md new file mode 100644 index 000000000..c2186350c --- /dev/null +++ b/examples/general_sampling/README.md @@ -0,0 +1,157 @@ +## Compilation + +Create a build directory. +Build the example by running the following commands in the build directory you created. + +```bash +cmake +make +``` +You might have to specify the path to liblpsolve55.so/dll/dylib. Try these: + +```bash +cmake . -DLP_SOLVE=_PATH_TO_LIB_FILE +make +``` +To find where this path is Try this command + + find /usr -name "liblpsolve55.so" 2>/dev/null + For example: -DLP_SOLVE=/usr/lib/lpsolve/liblpsolve55.so + +If you dont find it try downloading it. Most common way is: + + sudo apt update + sudo apt install lpsolve55 lpsolve55-dev + +## Usage: +```bash + ./general_sampling_time_limit +``` + +***How to use general_sampling_time_limit.cpp*** + +# Introduction + +The main focus of the example is to apply the uniform sampling methods of Volesti. +The polytope is generated by Volesti generators or is a custom Ax<=b polytope with A and b in csv form. +Then, most of the uniform methods are used to sample from it. The sampling lasts until an ESS limit +is reached or after you generate as many samples as you want. +In the end, the code returns (among others): + - Total samples + - ESS + - Time + +There are many custom parameters you can modify. Unfortunatelly, you currently need to modify the source code: + +# Choosing your desired sampling method + +You can choose to use one or all of the uniform sampling methods provided. +To do that check lines 592-598. Set true to the methods you want to run and false to the ones you don't. + +# Batch size: + +In line 138 you can find a function called compute_batch_size. There you can define how many samples each batch will return. +In order to reach a specific ess we sample in small batches and calculate the ESS. There are a few things to notice here: + - The smaller the batch size, the more time is spent on calculating the ESS. + - The higher the batch size the more difficult it is to reach the target exactly. + +Thus a compromise is needed. We provide numbers for some of the most common polytopes after many trial and error experiments. You can follow them or use your own. + +# Dynamic Batch size: + +A dynamic batch size calculation is also implemented like this: we begin with 10*target_ESS samples and check their ESS. +Calculating the efficiency of these samples, we request as many samplse as needed to reach the target ESS. After 3 batches at most target is reached. +You can turn on/off the dynamic batch size following these steps: + +How to turn ON dynamic batch size: + - Comment lines 271-272 and uncomment lines 273-274. + - Uncomment lines 337-352. +How to turn OFF dynamic batch size: + - Uncomment lines 271-272 and comment lines 273-274. + - Comment lines 337-352. + +Recommendation: Don't start your sampling with the dynamic batch size. There is a chance that the sampler will get stuck or that +it will take TOO long to sample 10 times your target. It is better (especially for an unknown polytope) to begin with a small +custom batch size, check the behaviour and then use the dynamic batch size is all is ok. + +# Walk lenght: + +In line 240 you cn find a set_walk_len function. Walk length works like a thinning factor. The values already set are standard. +You should not change then unless you want to experiment with different values if needed. Just change the return value to what you want. + +# Volume Shell Analysis: + +Uncomment lines 412-425 to see the full Shell Analysis of the KS test. These lines are commented by default to avoid clutter. + +# Print samples to file + +Uncomment lines 431-432 (same for lines 563-564 if you use CRHMC) to print output samples to a txt file for later use. +Carefull as only one file is kept, the last written, so if you aim for a specific method samples run only that method. + +# Choose native Polytope + +Volesti already provides generators for some classic polytopes. In lines 642-646 you can select one of them to sample from. +Just uncomment the one you want and comment the rest. You can select from Cube,Simplex,Cross,Birkhoff,Skinny Cube. + +# Custom Polytope + +If you want to sample from a custom polytope, we offer an option. The polytope must be in H-representation (A*x <= b). +A and b matrices must be in csv form. + +How to use custom mode: + - In line 611 you must set this boolean to true (false for native polytopes) + - In line 618 you must provide the names of the csv files holding A and b matrices. + +# Choose dimensions, target ESS and time limit + + - In line 582 there is a vector holding the dimensions. Add the dimensions you want the sampling to take place. We later loop this vector + so you can add multiple dimensions. (example {3,5,6,7,8,9,10}) + - In line 584 you can choose the target ESS. + - In order to sample only for as many samples as you choose, set target ESS to 1 and batch size to how many samples you want. + - In line 588 is the time limit in seconds. After some sampling methods exceed this limit, the method is aborted and the current results saved. + +# Rotate polytope + +We offer the chance to rotate any polytope around each plane. This is usefull to reduce bias from the fact that generating polytopes +sits aligned with the axis and thus helping some of the sampling methods. + +How to rotate a polytope + - In lines 648-649 choose an angle or define your own angle. + +# Print the polytope + +Sometimes, for debugging purposes it is useful to see the polytope that you sample from. + +How to print the polytope + - Uncomment or comment line 652 + +# Expected output (example) + +```text +Starting Benchmark. +Target ESS: 500 +Time Limit: 3600s per method. + +======================================== +*** Running for dimension 10 *** + +[AcceleratedBilliardWalk] Current batch_size: 2162 +[AcceleratedBilliardWalk] Current walk_len: 1 +[AcceleratedBilliardWalk] Samples: 10810 | ESS: 638 | Time: 0.248768s +[AcceleratedBilliardWalk] DONE. Final ESS: 638 +``` + +Results are also saved in a file after each sampling... +The file is called `benchmark_results.txt` and should look like this: + +```csv +Dim, Method, Time(s), Points, ESS +81, AcceleratedBilliardWalk, 0.248768, 10810, 638 +81, CDHRWalk, 0.188994, 6744, 662 +81, BilliardWalk, 1.12567, 40400, 570 +81, RDHRWalk, 26.7692, 59700, 552 +81, SparseBilliardWalk, 3.03094, 4324, 780 +``` + + + diff --git a/examples/general_sampling/general_sampling_time_limit.cpp b/examples/general_sampling/general_sampling_time_limit.cpp new file mode 100644 index 000000000..ec0f24320 --- /dev/null +++ b/examples/general_sampling/general_sampling_time_limit.cpp @@ -0,0 +1,795 @@ +#include +#include +#include +#include +#include +#include + +#include +#include "Eigen/Eigen" + +// These help with generating polytopes +#include "hpolytope.h" +#include "known_polytope_generators.h" +#include "custom_generators.h" +#include "order_polytope_generator.h" + +#include "cartesian_geom/cartesian_kernel.h" +#include "sampling/random_point_generators.hpp" +#include "random_walks/random_walks.hpp" +#include "random_walks/sparse_uniform_billiard_walk.hpp" +#include "preprocess/max_inscribed_ellipsoid.hpp" +#include "preprocess/inscribed_ellipsoid_rounding.hpp" +#include "convex_bodies/ellipsoid.h" +#include "convex_bodies/hpolytope.h" + +// These are related to PSRF and ESS and K-S test +#include "sampling/sample_correlation_matrices.hpp" +#include "matrix_operations/EigenvaluesProblems.h" +#include "diagnostics/effective_sample_size.hpp" +#include "diagnostics/univariate_psrf.hpp" +#include "diagnostics/scaling_ratio.hpp" +#include "diagnostics/KS_test.hpp" + +typedef double NT; +typedef Cartesian Kernel; +typedef typename Kernel::Point Point; +typedef Eigen::Matrix MT; +typedef Eigen::Matrix VT; +typedef BoostRandomNumberGenerator RNGType; +typedef HPolytope HPOLYTOPE; + +//Usefull struct to rotate polytope +template +HPOLYTOPE rotate_all_dims(const HPOLYTOPE& P, typename HPOLYTOPE::NT angle) +{ + using NT = typename HPOLYTOPE::NT; + int dim = P.dimension(); + + // Build global rotation matrix + Eigen::Matrix R = + Eigen::Matrix::Identity(dim, dim); + + NT c = std::cos(angle); + NT s = std::sin(angle); + + // Apply rotation in each adjacent coordinate plane + for (int k = 0; k < dim - 1; ++k) { + Eigen::Matrix Rk = + Eigen::Matrix::Identity(dim, dim); + + Rk(k, k) = c; + Rk(k, k+1) = -s; + Rk(k+1, k) = s; + Rk(k+1, k+1) = c; + + R = R * Rk; // compose rotations + } + + // Extract A and b + auto A = P.get_mat(); + auto b = P.get_vec(); + + // Apply A' = A R + Eigen::Matrix A_rot = A * R; + + return HPOLYTOPE(dim, A_rot, b); +} + + +//Used to print the ineqialities (if needed) +void print_hpoly(const HPOLYTOPE& P) { + const auto& A = P.get_mat(); // Matrix: rows = constraints, cols = dimension + const auto& b = P.get_vec(); // Vector: one entry per inequality + + unsigned int m = A.rows(); + unsigned int d = A.cols(); + + std::cout << "H-representation: A x <= b\n"; + std::cout << "Number of inequalities: " << m << "\n"; + std::cout << "Dimension: " << d << "\n\n"; + + for (unsigned int i = 0; i < m; ++i) { + std::cout << i << ": "; + for (unsigned int j = 0; j < d; ++j) { + std::cout << A(i, j) << " * x" << j; + if (j < d - 1) std::cout << " + "; + } + std::cout << " <= " << b(i) << "\n"; + } + std::cout << std::endl; +} + + +// Used to write final samples to a file +void write_to_file(std::string filename, std::vector const& randPoints) { + std::ofstream out(filename); + auto coutbuf = std::cout.rdbuf(out.rdbuf()); //save and redirect + for(int i=0; i(end_time - start_time).count(); + total_time += elapsed; + if (!label.empty()) + std::cout << "[" << walk_name << "] " << label << " = " << elapsed << " s\n"; + return elapsed; + } + + double get_total_time() const { return total_time; } + +private: + std::string walk_name; + std::chrono::steady_clock::time_point start_time; + double total_time; +}; + +// Function to dynamically compute the batch size based on the walk type and dimension +unsigned int compute_batch_size(const std::string& walk_name, unsigned int dim) { + + // Default base size + unsigned int batch_size = 1000; + + // Adaptive logic by walk type + if (walk_name == "BallWalk") { + //Cube + //batch_size = dim * 150 + 1500; + //Simplex + //batch_size = dim * 500 + 1500; + //Birkhoff + batch_size = dim * 1000 + 1500; + } + else if (walk_name == "AcceleratedBilliardWalk") { + //Cube + //batch_size = dim * 2 + 200; + //Simplex + //batch_size = dim * 1 + 1500; + //Birkhoff + batch_size = dim * 4 + 2000; + } + else if (walk_name == "BilliardWalk") { + //Cube + //batch_size = dim * 4 + 300; + //Simplex + //batch_size = dim * 2 + 3000; + //Birkhoff + batch_size = dim * 20 + 2500; + } + else if (walk_name == "SparseBilliardWalk") { + //Cube + //batch_size = dim * 2 + 200; + //Simplex + //batch_size = dim * 1 + 1500; + //Birkhoff + batch_size = dim * 2 + 2000; + } + else if (walk_name == "CDHRWalk") { + //Cube + //batch_size = dim * 1 + 500; + //Simplex + //batch_size = dim * 2 + 2400; + //Birkhoff + batch_size = dim * 3 + 2400; + //Biology + //batch_size = dim * 100 + 2400; + } + else if (walk_name == "RDHRWalk") { + //Cube + //batch_size = dim * 15 + 500; + //Simplex + //batch_size = dim * 350 + 1500; + //Birkhoff + batch_size = dim * 350 + 1500; + } + else if (walk_name == "DikinWalk") { + batch_size = dim * 100 + 1000; + } + else if (walk_name == "JohnWalk") { + batch_size = dim * 250 + 1000; + } + else if (walk_name == "VaidyaWalk") { + batch_size = dim * 150 + 1000; + } + + + return batch_size; +} + +// Converts a vector of Points into an Eigen matrix +template +MT vector_to_eigen(const std::vector& someSamples) { + MT samples(someSamples[0].dimension(), someSamples.size()); + for (unsigned int jj = 0; jj < someSamples.size(); ++jj) + samples.col(jj) = someSamples[jj].getCoefficients(); + return samples; +} + +// Computes ESS on a given Eigen matrix and returns the min ESS +template +unsigned int compute_ess(const MT& samples) { + unsigned int min_ess = 0; + VT ess_vector = effective_sample_size(samples, min_ess); + //std::cout << "Current ESS is: " << min_ess << "\n"; + return min_ess; +} + +// Computes PSRF for all accumulated points (no timing inside) +template +double compute_psrf(const std::vector& someSamples) { + // Convert to Eigen matrix + MT finalSamples(someSamples[0].dimension(), someSamples.size()); + for (unsigned int jj = 0; jj < someSamples.size(); ++jj) + finalSamples.col(jj) = someSamples[jj].getCoefficients(); + + // Compute PSRF + VT psrf = univariate_psrf(finalSamples); + double max_psrf = psrf.maxCoeff(); + + return max_psrf; +} + +// We choose the walk_len appropriate for each method +unsigned int set_walk_len(const std::string& walk_name, unsigned int dim) { + static const std::unordered_map> walk_len_map = { + {"BallWalk", [](unsigned int dim){ return dim*2; }}, + {"BilliardWalk", [](unsigned int){ return 1; }}, + {"AcceleratedBilliardWalk", [](unsigned int){ return 1; }}, + {"CDHRWalk", [](unsigned int dim){ return dim*2; }}, + {"RDHRWalk", [](unsigned int dim){ return dim*2; }}, + {"DikinWalk", [](unsigned int dim){ return 50+dim; }}, + {"JohnWalk", [](unsigned int dim){ return 50+dim; }}, + {"VaidyaWalk", [](unsigned int dim){ return 50+dim; }} + }; + + if (auto it = walk_len_map.find(walk_name); it != walk_len_map.end()) + return it->second(dim); + + return 1; // default fallback +} + +template +bool sample_using_walk(HPOLYTOPE& Polytope, + Point const& start_point, + RNGType& rng, + unsigned int target_ESS, + const std::string& walk_name, + double time_limit_sec) +{ + Point starting_point = start_point; + unsigned int dim = Polytope.dimension(); + + Timer t; + //Uncommenct these 2 lines and comment lines 269-270 and 334-349 to swap back to manual batch size + unsigned int batch_size = compute_batch_size(walk_name, dim); + std::cout << "\n[" << walk_name << "] Current batch_size: " << batch_size << "\n"; + //unsigned int batch_size = target_ESS*5; + //std::cout << "\n[" << walk_name << "] Using dynamic batch size. Initial batch_size: " << batch_size << "\n"; + + unsigned int walk_len = set_walk_len(walk_name, dim); + std::cout << "[" << walk_name << "] Current walk_len: " << walk_len << "\n"; + + unsigned int current_ESS = 0; + unsigned int loop_step = 1; + bool failed_to_converge = false; + bool timed_out = false; + + Timer eigen_timer(walk_name); + Timer ess_timer(walk_name); + Timer generator_timer(walk_name); + Timer psrf_timer(walk_name); + + typedef RandomPointGenerator Generator; + std::vector allSamples; + allSamples.reserve(target_ESS * 10); + + while (current_ESS < target_ESS) { + + // CHECK TIME LIMIT + if (generator_timer.get_total_time() > time_limit_sec) { + std::cout << "[" << walk_name << "] TIMEOUT (" + << generator_timer.get_total_time() << "s > " + << time_limit_sec << "s). Stopping." + << std::string(20, ' ') << "\n"; // Adds 20 spaces to clear the line; + timed_out = true; + failed_to_converge = true; + break; + } + + //unsigned int num_points = batch_size; + std::vector batchPoints; + + generator_timer.start(); + Generator::apply(Polytope, starting_point, batch_size, walk_len, batchPoints, push_back_policy, rng); + generator_timer.stop(""); // Accumulates time + + if (!batchPoints.empty()) { + starting_point = batchPoints.back(); + } + allSamples.insert(allSamples.end(), batchPoints.begin(), batchPoints.end()); + + // Check ESS + if (allSamples.size() >= target_ESS) { + eigen_timer.start(); + MT samples = vector_to_eigen(allSamples); + eigen_timer.stop(""); + + ess_timer.start(); + current_ESS = compute_ess(samples); + ess_timer.stop(""); + + if (current_ESS >= target_ESS) break; + + // Optional: Print progress + std::cout << "[" << walk_name << "] Samples: " << allSamples.size() + << " | ESS: " << current_ESS + << " | Time: " << generator_timer.get_total_time() << "s\r" << std::flush; + } else { + current_ESS = 0; + } + + // // // THE DYNAMIC batch size logic + // // Calculate how useful each sample is + // double ess_per_sample = (double)current_ESS / (double)allSamples.size(); + + // if (ess_per_sample > 1e-6) { + // unsigned int remaining_ESS = target_ESS - current_ESS; + // // Predict needed samples + 10% safety buffer + // unsigned int samples_needed = static_cast((remaining_ESS / ess_per_sample) * 1.1); + + // // Don't let the batch size get smaller than the initial guess, + // // but allow it to scale up to finish the job in the next loop. + // batch_size = std::max(batch_size, samples_needed); + // } else { + // // If efficiency is essentially zero, double the batch to try and find signal + // batch_size *= 2; + // } + + // Failsafe for infinite loops + if (loop_step > 50 && current_ESS < target_ESS) { + std::cout << "\n[" << walk_name << "] I am sorry. I did not converge!\n"; + failed_to_converge = true; + break; + } + loop_step++; + } + std::cout << "[" << walk_name << "] Samples: " << allSamples.size() + << " | ESS: " << current_ESS + << " | Time: " << generator_timer.get_total_time() << "s\r" << std::flush; + + // PRINT RESULTS TO CONSOLE + if(!failed_to_converge) { + std::cout << "\n[" << walk_name << "] DONE. Final ESS: " << current_ESS << "\n"; + std::cout << "[" << walk_name << "] Total generation time: " << generator_timer.get_total_time() << " s\n"; + } + + // SAVE TO FILE + // Format: Dimension, WalkName, Time, Points, ESS + std::ofstream outfile; + outfile.open("benchmark_results.txt", std::ios_base::app); // Append mode + + if (outfile.is_open()) { + outfile << Polytope.dimension() << ", " + << walk_name << ", " + << generator_timer.get_total_time() << ", " + << allSamples.size() << ", " + << current_ESS << "\n"; + outfile.close(); + std::cout << "[" << walk_name << "] Results saved to benchmark_results.txt" + << std::string(20, ' ') << "\n"; // Adds 20 spaces to clear the line; + } else { + std::cerr << "Unable to open file to write results!\n"; + } + + // PSRF calculation + psrf_timer.start(); + double max_psrf = compute_psrf(allSamples); + psrf_timer.stop("Total PSRF time"); + + // PSRF result + std::cout << "[" << walk_name << "] Total time to calculate ESS = " << ess_timer.get_total_time() << " s\n"; + std::cout << "[" << walk_name << "] PSRF = " << max_psrf << std::endl; + + // K-S statistical test //////////////////////////////////////////////////// + MT samples_mat = vector_to_eigen(allSamples); + int computed_thin = static_cast(samples_mat.cols() / current_ESS); + int thin_factor = std::max(10, computed_thin * 2); + if (thin_factor < 1) thin_factor = 1; + + auto [ks_stat, p_val, observed, expected] = global_scaling_test(Polytope, samples_mat, thin_factor); + + std::cout << "[" << walk_name << "] KS Statistic: " << ks_stat << "\n"; + std::cout << "[" << walk_name << "] P-Value: " << p_val << "\n"; + + // Comment-uncomment to see full volume analysis + // Print the Shell Analysis Table (using the unpacked variables) + // std::cout << "Volume Shells Analysis (Expected vs Observed):\n"; + // std::cout << "Exp Vol% | Obs Vol% | Deviation\n"; + // std::cout << "-------------------------------\n"; + + // double max_dev_perc = 0.0; + // for(size_t i = 0; i < 10; ++i) { + // double dev = (observed[i] - expected[i]) * 100.0; + + // if(std::abs(dev) > std::abs(max_dev_perc)) max_dev_perc = dev; + + // printf(" %4.1f%% | %4.1f%% | %+6.2f%%\n", + // expected[i]*100.0, observed[i]*100.0, dev); + // } + // std::cout << "-------------------------------\n"; + + if (timed_out || failed_to_converge) return false; + + // Uncomment to write samples to file + //write_to_file("All_Samples.txt", allSamples); + //std::cout << "All samples written to ALL_Samples.txt file. DONE" << std::endl; + + return true; +} + +///////////////RIEMMANIAN////////////////////////////////////////////////////////////////////////////////////////////////////// +template +bool sample_using_crhmc(PolytopeType& HP, + typename PolytopeType::PointType& /*center*/, + RNGType& rng, + unsigned int target_ESS, + const std::string& walk_name) +{ + using NT = double; + using Point = typename PolytopeType::PointType; + using VT = Eigen::Matrix; + using MT = typename PolytopeType::MT; + + using Func = ZeroScalarFunctor; + using Grad = ZeroFunctor; + using Hess = ZeroFunctor; + + // --- Configuration --- + double current_efficiency = 1.0 / 10.0; + int n_burns = 1000; + int walk_len = 1; // Thinning + + // Global container for ALL samples across batches + std::vector all_samples; + all_samples.reserve(target_ESS * 20); + + Timer total_timer(walk_name); + total_timer.start(); + + double current_ESS = 0.0; + int batch_count = 1; + + std::cout << "[" << walk_name << "] Target ESS: " << target_ESS << "\n"; + + // --- The Smart Loop --- + while (current_ESS < target_ESS) { + + // 1. Calculate how many samples we need + double missing_ESS = target_ESS - current_ESS; + + // "Smart Batching": Estimate samples needed based on current efficiency + // We add a 10% buffer (1.1) to try and finish in this batch + int n_samples_needed = static_cast((missing_ESS / current_efficiency) * 1.1); + + // Safety clamps + if (n_samples_needed < 1000) n_samples_needed = 1000; + if (n_samples_needed > 100000) n_samples_needed = 100000; // Cap to prevent memory explosion + + std::cout << "[" << walk_name << "][Batch " << batch_count << "] Requesting " + << n_samples_needed << " samples (Efficiency: " << current_efficiency << ")\n"; + + // Setup Helper Objects (Must be fresh per run) + Func* f = new Func; + Grad* g = new Grad; + std::list batch_list; + + // EXECUTE CRHMC + // Note: passing 'rng' ensures the random sequence continues, making this valid + execute_crhmc, Grad, Func, Hess, CRHMCWalk, 1>( + HP, rng, batch_list, walk_len, n_samples_needed, n_burns, g, f + ); + + delete f; + delete g; + + // Merge Samples + // We move elements from list to vector to avoid copying + all_samples.insert(all_samples.end(), batch_list.begin(), batch_list.end()); + + // 5. Update Statistics + // We must convert ALL samples to matrix to calculate total ESS + MT samples_matrix = MT(HP.dimension(), all_samples.size()); + for (size_t i = 0; i < all_samples.size(); ++i) { + samples_matrix.col(i) = all_samples[i].getCoefficients(); + } + + // Calculate ESS + current_ESS = compute_ess(samples_matrix); + + // Update Efficiency for next loop + // Efficiency = ESS / Total_Raw_Samples + if (all_samples.size() > 0) { + current_efficiency = current_ESS / all_samples.size(); + } + + std::cout << "[" << walk_name << "][Batch " << batch_count << "] Current ESS: " << current_ESS + << " / " << target_ESS << "\n"; + + // Break if we are stuck (Efficiency drops too low) + if (current_efficiency < 0.0001 && all_samples.size() > 5000) { + std::cout << "[" << walk_name << "] CRITICAL: Efficiency too low. Stopping.\n"; + break; + } + + batch_count++; + } + + total_timer.stop(""); + + // --- Final Reporting & Tests --- + + // 1. PSRF + Timer psrf_timer(walk_name); + psrf_timer.start(); + double max_psrf = compute_psrf(all_samples); + psrf_timer.stop(""); + + // 2. KS Test + MT final_matrix = MT(HP.dimension(), all_samples.size()); + for (size_t i = 0; i < all_samples.size(); ++i) { + final_matrix.col(i) = all_samples[i].getCoefficients(); + } + + // KS Logic + double safe_ess = (current_ESS > 0.0) ? current_ESS : 1.0; + int computed_thin = static_cast(final_matrix.cols() / safe_ess); + int thin_factor = std::max(10, computed_thin * 2); + if (thin_factor < 1) thin_factor = 1; + + auto [ks_stat, p_val, observed, expected] = global_scaling_test(HP, final_matrix, thin_factor); + + std::cout << "------------------------------------------------\n"; + std::cout << "[" << walk_name << "] Total time : " << total_timer.get_total_time() << " s\n"; + std::cout << "[" << walk_name << "] Total Batches : " << (batch_count - 1) << "\n"; + std::cout << "[" << walk_name << "] Total Samples : " << all_samples.size() << "\n"; + std::cout << "[" << walk_name << "] Final ESS : " << current_ESS << "\n"; + std::cout << "[" << walk_name << "] Final PSRF : " << max_psrf << "\n"; + std::cout << "[" << walk_name << "] KS Statistic : " << ks_stat << "\n"; + std::cout << "[" << walk_name << "] P-Value : " << p_val << "\n"; + std::cout << "------------------------------------------------\n"; + + // Uncomment to write samples to file + //write_to_file("All_Samples.txt", all_samples); + //std::cout << "All samples written to ALL_Samples.txt file. DONE" << std::endl; + return true; +} +/////////////////End of Riemannian/////////////////////////////////////////////////////////////////////////////////////////////////////// + +typedef BallWalk::template Walk BallWalkType; +typedef BilliardWalk::template Walk BilliardWalkType; +typedef AcceleratedBilliardWalk::template Walk AcceleratedBilliardWalkType; +typedef CDHRWalk::template Walk CDHRWalkType; +typedef DikinWalk::template Walk DikinWalkType; +typedef JohnWalk::template Walk JohnWalkType; +typedef RDHRWalk::template Walk RDHRWalkType; +typedef VaidyaWalk::template Walk VaidyaWalkType; +typedef SparseBilliardWalk::template Walk SparseBilliardaWalkType; + +int main(int argc, char const *argv[]) { + + // You can adjust dimensions as needed + std::vector dimensions = {25}; + // Select target ESS + unsigned int target_ESS = 800; + + //seed for order polytopes + int base_seed = 42; + + // Time Limit: 1 Hour (3600 seconds) + // If a method exceeds this, it stops and is skipped for all larger dimensions. + double TIME_LIMIT_SEC = 1200.0; + + // Choose what methods to use. True is used, false is not used. + std::map active_methods; + active_methods["AcceleratedBilliardWalk"] = true; + active_methods["CDHRWalk"] = false; + active_methods["BallWalk"] = false; + active_methods["BilliardWalk"] = false; + active_methods["RDHRWalk"] = false; + active_methods["CRHMCWalk"] = false; + active_methods["SparseBilliardWalk"] = false; + + // Initialize output file (overwrite old one) + std::ofstream outfile("benchmark_results.txt"); + outfile << "Dim, Method, Time(s), Points, ESS\n"; + outfile.close(); + + std::cout << "Starting Benchmark.\n"; + std::cout << "Target ESS: " << target_ESS << "\n"; + std::cout << "Time Limit: " << TIME_LIMIT_SEC << "s per method.\n"; + + //////CUSTOM Polytopes in A*x<=b form/////////////*******************************************//////////////////////////////////////// + // TOGGLE THIS: Set to 'true' for your custom CSVs, 'false' for the Cube/Simplex/... benchmark + bool USE_CUSTOM_MODEL = true; + + HPOLYTOPE custom_polytope; // Placeholder for the loaded model + + if (USE_CUSTOM_MODEL) { + try { + // Load the model ONCE before the loop. Place the csv files in the build folder. + custom_polytope = load_custom_polytope("Birkhoff25_volumetric_barrier_A.csv", "Birkhoff25_volumetric_barrier_b.csv"); + + // Overwrite dimensions list to run exactly ONCE for the model's dimension + dimensions = { static_cast(custom_polytope.dimension()) }; + + std::cout << ">>> CUSTOM MODE ACTIVATED: Loaded model with " << dimensions[0] << " dimensions.\n"; + } catch (const std::exception& e) { + std::cerr << "CRITICAL ERROR: Could not load custom files. " << e.what() << std::endl; + return 1; + } + } + //////End of Custom Polytopes///////////////////*******************************************//////////////////////////////////////// + + for (auto dim : dimensions) { + + std::cout << "\n" << std::string(40, '=') << "\n"; + std::cout << "*** Running for dimension " << dim << " ***\n"; + + HPOLYTOPE Polytope_simple; + HPOLYTOPE Polytope; + if (USE_CUSTOM_MODEL) { + Polytope = custom_polytope; + } else { + // Generate desired Polytope + //Polytope_simple = generate_cube(dim, false); + Polytope_simple = generate_birkhoff(dim); + //Polytope_simple = generate_cross(dim, false); + //Polytope_simple = generate_skinny_cube(dim,false); + //Polytope_simple = generate_simplex(dim, false); + + // Generate order polytopes + // unsigned int m = 3 * dim; + // int current_seed = base_seed + dim; + // std::cout << "\nCreating order polytope...\n"; + // Polytope_simple = random_orderpoly(dim, m, current_seed); + + //double angle = 53.0 * M_PI / 180.0; //53 deg + double angle = 0.0 * M_PI / 180.0; //0 deg + Polytope = rotate_all_dims(Polytope_simple, angle); + + //print_hpoly(Polytope); //Uncomment to print the polytope + } + // Setup RNG and Starting Point + RNGType rng(Polytope.dimension()); + auto inner = Polytope.ComputeInnerBall(); + Point center = inner.first; + + //Use the following lines to force the center (starting point) be 0 + // point> center(dim); + // center.set_to_origin(); + // std::cout << "Starting point: "; + // center.print(); + + // ========================================================= + // ACCELERATED BILLIARD WALK + // ========================================================= + if (active_methods["AcceleratedBilliardWalk"]) { + bool success = sample_using_walk( + Polytope, center, rng, target_ESS, "AcceleratedBilliardWalk", TIME_LIMIT_SEC + ); + if (!success) { + active_methods["AcceleratedBilliardWalk"] = false; + std::cout << "!!! Disabling AcceleratedBilliardWalk for future dimensions.\n"; + } + } else { + std::cout << "\n[AcceleratedBilliardWalk] Skipping (previously timed out).\n"; + } + + // ========================================================= + // CDHR (Coordinate Directions Hit-and-Run) + // ========================================================= + if (active_methods["CDHRWalk"]) { + bool success = sample_using_walk( + Polytope, center, rng, target_ESS, "CDHRWalk", TIME_LIMIT_SEC + ); + if (!success) { + active_methods["CDHRWalk"] = false; + std::cout << "!!! Disabling CDHR for future dimensions.\n"; + } + } else { + std::cout << "\n[CDHRWalk] Skipping (previously timed out).\n"; + } + + // ========================================================= + // BILLIARD WALK + // ========================================================= + if (active_methods["BilliardWalk"]) { + bool success = sample_using_walk( + Polytope, center, rng, target_ESS, "BilliardWalk", TIME_LIMIT_SEC + ); + if (!success) { + active_methods["BilliardWalk"] = false; + std::cout << "!!! Disabling Billiard Walk for future dimensions.\n"; + } + } else { + std::cout << "\n[BilliardWalk] Skipping (previously timed out).\n"; + } + + // ========================================================= + // RDHR (Random Directions Hit-and-Run) + // ========================================================= + if (active_methods["RDHRWalk"]) { + bool success = sample_using_walk( + Polytope, center, rng, target_ESS, "RDHRWalk", TIME_LIMIT_SEC + ); + if (!success) { + active_methods["RDHRWalk"] = false; + std::cout << "!!! Disabling RDHR for future dimensions.\n"; + } + } else { + std::cout << "\n[RDHRWalk] Skipping (previously timed out).\n"; + } + + // ========================================================= + // BALL WALK + // ========================================================= + if (active_methods["BallWalk"]) { + bool success = sample_using_walk( + Polytope, center, rng, target_ESS, "BallWalk", TIME_LIMIT_SEC + ); + if (!success) { + active_methods["BallWalk"] = false; + std::cout << "!!! Disabling Ball Walk for future dimensions.\n"; + } + } else { + std::cout << "\n[BallWalk] Skipping (previously timed out).\n"; + } + + // ========================================================= + // SparseBilliardWalk + // ========================================================= + if (active_methods["SparseBilliardWalk"]) { + bool success = sample_using_walk( + Polytope, center, rng, target_ESS, "SparseBilliardWalk", TIME_LIMIT_SEC + ); + if (!success) { + active_methods["SparseBilliardWalk"] = false; + std::cout << "!!! Disabling Sparse Billiard Walk for future dimensions.\n"; + } + } else { + std::cout << "\n[SparseBilliardWalk] Skipping (previously timed out).\n"; + } + + // ========================================================= + // CRHMC WALK + // ========================================================= + if (active_methods["CRHMCWalk"]) { + std::cout << "Starting CRHMC..." << std::endl; + bool success = sample_using_crhmc(Polytope, center, rng, target_ESS, "CRHMCWalk"); + if (!success) { + active_methods["CRHMCWalk"] = false; + std::cout << "!!! Disabling CRHMCWalk for future dimensions.\n"; + } + } else { + std::cout << "\n[CRHMCWalk] Skipping (previously timed out).\n"; + } + + } // End dimension loop + + std::cout << "\nBenchmark Complete.\n"; + return 0; +} \ No newline at end of file diff --git a/include/diagnostics/KS_test.hpp b/include/diagnostics/KS_test.hpp new file mode 100644 index 000000000..af89311e6 --- /dev/null +++ b/include/diagnostics/KS_test.hpp @@ -0,0 +1,134 @@ +// VolEsti (volume computation and sampling library) + +//KS Test (radial distribution) + +#ifndef DIAGNOSTICS_KS_TEST_HPP +#define DIAGNOSTICS_KS_TEST_HPP + +#include +#include +#include +#include +#include + +// Kolmogorov distribution: P(K > z) +inline double kolmogorov_prob(double z) { + if (z <= 0.0) return 1.0; + double sum = 0.0; + for (int k = 1; k <= 200; ++k) { + double term = std::exp(-2.0 * k * k * z * z); + sum += (k % 2 ? term : -term); + if (term < 1e-15) break; + } + return std::max(0.0, std::min(1.0, 2.0 * sum)); +} + +// Global uniformity test +template +std::tuple, std::vector> +global_scaling_test(const Polytope& P, + const typename Polytope::MT& samples, + int thinning_factor = 1) +{ + using VT = typename Polytope::VT; + + const int dim = P.dimension(); + const int n_total = static_cast(samples.cols()); + + // Setup center and constraints + VT center = samples.rowwise().mean(); + const auto A = P.get_mat(); + const auto b = P.get_vec(); + VT b_shifted = b - A * center; + + if (b_shifted.minCoeff() < 1e-12) { + std::cerr << "[GlobalKS] Warning: Empirical center is on boundary/outside.\n"; + } + std::vector rvals_all; + rvals_all.reserve(n_total); + + for (int i = 0; i < n_total; ++i) { + VT q = samples.col(i) - center; + VT u = A * q; + + double r_max = 0.0; + for (int k = 0; k < u.size(); ++k) { + if (u[k] > 0.0) { + double denom = b_shifted[k]; + if (denom > 1e-14) { + double t = u[k] / denom; + if (t > r_max) r_max = t; + } else { + r_max = 1.0; + } + } + } + if (r_max > 1.0) r_max = 1.0; + if (r_max < 0.0) r_max = 0.0; + rvals_all.push_back(r_max); + } + + // We only perform the KS test on the subset to ensure independence. + std::vector uvals; + if (thinning_factor < 1) thinning_factor = 1; + + // Reserve roughly N / factor + uvals.reserve(n_total / thinning_factor + 1); + + for (int i = 0; i < n_total; i += thinning_factor) { + double r = rvals_all[i]; + // Transform r -> u = r^d (Probability Integral Transform) + uvals.push_back(std::pow(r, dim)); + } + + // Compute KS Statistic + std::sort(uvals.begin(), uvals.end()); + int N_test = static_cast(uvals.size()); + + if (N_test < 5) return {0.0, 1.0, {}, {}}; // Too few samples + + double ks_stat = 0.0; + for (int i = 0; i < N_test; ++i) { + double F_emp_lo = static_cast(i) / N_test; + double F_emp_hi = static_cast(i + 1) / N_test; + double F_theo = uvals[i]; + + double diff = std::max(std::abs(F_emp_lo - F_theo), + std::abs(F_emp_hi - F_theo)); + if (diff > ks_stat) ks_stat = diff; + } + + // Standard P-value + double sqrt_n = std::sqrt(static_cast(N_test)); + double lambda = (sqrt_n + 0.12 + 0.11 / sqrt_n) * ks_stat; + double p_value = kolmogorov_prob(lambda); + + // Shell diagnostics + std::vector exp_coverage(10), obs_coverage(10, 0.0); + std::vector r_thresholds(10); + std::vector shell_counts(10, 0); + + for (int k = 0; k < 10; ++k) { + exp_coverage[k] = 0.1 * (k + 1); + r_thresholds[k] = std::pow(exp_coverage[k], 1.0 / dim); + } + + for (double r : rvals_all) { + for (int k = 0; k < 10; ++k) { + if (r <= r_thresholds[k]) { + shell_counts[k]++; + break; + } + } + } + + int cumulative = 0; + for (int k = 0; k < 10; ++k) { + cumulative += shell_counts[k]; + obs_coverage[k] = static_cast(cumulative) / n_total; + } + + return {ks_stat, p_value, obs_coverage, exp_coverage}; +} + +#endif diff --git a/include/generators/custom_generators.h b/include/generators/custom_generators.h new file mode 100644 index 000000000..13ce24b1c --- /dev/null +++ b/include/generators/custom_generators.h @@ -0,0 +1,83 @@ +#ifndef CUSTOM_GENERATORS_HPP +#define CUSTOM_GENERATORS_HPP + +#include +#include +#include +#include +#include +#include + +// ------------------------------------------------------------------------- +// Helper: Reads a CSV file into an Eigen Matrix +// ------------------------------------------------------------------------- +template +Eigen::Matrix read_csv_to_eigen(const std::string &path) { + std::ifstream indata; + indata.open(path); + + if (!indata.is_open()) { + throw std::runtime_error("Could not open file: " + path); + } + + std::string line; + std::vector values; + unsigned int rows = 0; + + while (std::getline(indata, line)) { + std::stringstream lineStream(line); + std::string cell; + while (std::getline(lineStream, cell, ',')) { + // Check for empty cells usually caused by trailing commas + if (!cell.empty()) { + values.push_back(static_cast(std::stod(cell))); + } + } + ++rows; + } + + if (rows == 0) return Eigen::Matrix(); + + // Calculate columns + unsigned int cols = values.size() / rows; + + // Map the std::vector to an Eigen Matrix + // We use RowMajor because CSVs are read row by row + return Eigen::Map>(values.data(), rows, cols); +} + +// ------------------------------------------------------------------------- +// Generator: Loads a Polytope defined by Ax <= b from CSV files +// ------------------------------------------------------------------------- +template +Polytope load_custom_polytope(const std::string &file_A, const std::string &file_b) { + + // Define types based on the Polytope template + typedef typename Polytope::NT NT; // Number Type (e.g., double) + typedef typename Eigen::Matrix MT; // Matrix Type + typedef typename Polytope::VT VT; // Vector Type + + std::cout << "Loading polytope from CSVs..." << std::endl; + + // Load matrices using the helper + MT A_raw = read_csv_to_eigen(file_A); + MT b_raw = read_csv_to_eigen(file_b); + + // Safety check + if (A_raw.rows() != b_raw.rows()) { + throw std::runtime_error("Dimension mismatch: Rows in A do not match rows in b."); + } + + // Convert b from Matrix (Nx1) to Vector (N) + VT b = b_raw.col(0); + + unsigned int dim = A_raw.cols(); + unsigned int num_constraints = A_raw.rows(); + + std::cout << "Successfully loaded: " << num_constraints << " constraints in " << dim << " dimensions." << std::endl; + + // Return the Polytope + return Polytope(dim, A_raw, b); +} + +#endif // CUSTOM_GENERATORS_HPP \ No newline at end of file diff --git a/include/random_walks/boundary_cdhr_walk.hpp b/include/random_walks/boundary_cdhr_walk.hpp index 5e6b5224c..31245717a 100644 --- a/include/random_walks/boundary_cdhr_walk.hpp +++ b/include/random_walks/boundary_cdhr_walk.hpp @@ -62,6 +62,8 @@ struct BCDHRWalk p2.set_coord(_rand_coord, _p_prev[_rand_coord] + bpair.second); } + const Point& getCurrentPoint() const noexcept { return _p; } + private : template diff --git a/include/random_walks/boundary_rdhr_walk.hpp b/include/random_walks/boundary_rdhr_walk.hpp index 29d500b6c..6fe843c19 100644 --- a/include/random_walks/boundary_rdhr_walk.hpp +++ b/include/random_walks/boundary_rdhr_walk.hpp @@ -57,6 +57,8 @@ struct BRDHRWalk } } + const Point& getCurrentPoint() const noexcept { return _p; } + private : template diff --git a/include/random_walks/random_walks.hpp b/include/random_walks/random_walks.hpp index 27fb39d6a..4c8cb5794 100644 --- a/include/random_walks/random_walks.hpp +++ b/include/random_walks/random_walks.hpp @@ -32,4 +32,6 @@ #include "random_walks/nuts_hmc_walk.hpp" #include "random_walks/langevin_walk.hpp" #include "random_walks/crhmc/crhmc_walk.hpp" +#include "random_walks/shake_and_bake_walk.hpp" +#include "random_walks/billiard_shake_and_bake_walk.hpp" #endif // RANDOM_WALKS_RANDOM_WALKS_HPP diff --git a/include/random_walks/sparse_uniform_billiard_walk.hpp b/include/random_walks/sparse_uniform_billiard_walk.hpp index 9b9df8bdb..bdf3be438 100644 --- a/include/random_walks/sparse_uniform_billiard_walk.hpp +++ b/include/random_walks/sparse_uniform_billiard_walk.hpp @@ -18,6 +18,7 @@ #include "convex_bodies/hpolytope.h" #include "sampling/sphere.hpp" #include "generators/boost_random_number_generator.hpp" +#include "preprocess/barrier_center_ellipsoid.hpp" struct SparseBilliardWalk { @@ -54,6 +55,52 @@ struct Walk typedef Eigen::Matrix MT; SparseRowMT _A_original; + + // Custom Constructor + template + Walk(GenericPolytope& P, const Point& p, RandomNumberGenerator& rng) + { + using MT = Eigen::Matrix; + using VT = Eigen::Matrix; + + _Len = NT(6.0) * std::sqrt(static_cast(P.dimension())); + + MT A_dense = MT(P.get_mat()); + VT b_dense = P.get_vec(); + + auto result = barrier_center_ellipsoid_linear_ineq( + A_dense, b_dense + ); + + MT H_dense = std::get<0>(result); + bool converged = std::get<2>(result); + + if (!converged) { + VT p_coeffs = p.getCoefficients(); + VT slack = b_dense - A_dense * p_coeffs; + MT S_inv_sq = MT::Zero(P.num_of_hyperplanes(), P.num_of_hyperplanes()); + for (unsigned int i = 0; i < P.num_of_hyperplanes(); ++i) { + NT s_val = (slack(i) < 1e-12) ? 1e-12 : slack(i); + S_inv_sq(i, i) = 1.0 / (s_val * s_val); + } + H_dense = A_dense.transpose() * S_inv_sq * A_dense; + } + + SparseMT H = H_dense.sparseView(); + + compute_cholesky_and_transformations(H); + + _b = P.get_vec(); + _A_original = P.get_mat().sparseView(); + _oracle_params.emplace(_L_inv, _A_original, _b); + + VT p_original = p.getCoefficients(); + VT p_rounded = _L_inv.transpose().template triangularView() * p_original; + Point p_rounded_point(p_rounded); + + initialize(P, p_rounded_point, rng); + } + // End of custom constructor template Walk(GenericPolytope& P, const Point& p, RandomNumberGenerator& rng, @@ -258,4 +305,4 @@ struct Walk }; }; -#endif // RANDOM_WALKS_SPARSE_BILLIARD_WALK_HPP \ No newline at end of file +#endif // RANDOM_WALKS_SPARSE_BILLIARD_WALK_HPP diff --git a/include/sampling/sampling.hpp b/include/sampling/sampling.hpp index e25494d1b..1603212c3 100644 --- a/include/sampling/sampling.hpp +++ b/include/sampling/sampling.hpp @@ -405,9 +405,9 @@ void crhmc_sampling(PointList &randPoints, NT, NegativeGradientFunctor > walk_params; - Point p = Point(problem.center); + Point p = Point(problem.center); problem.options.simdLen=simdLen; - walk_params params(input.df, p.dimension(), problem.options); + walk_params params(input.df, p.dimension(), problem.options); if (input.df.params.eta > 0) { params.eta = input.df.params.eta; @@ -415,15 +415,15 @@ void crhmc_sampling(PointList &randPoints, PushBackWalkPolicy push_back_policy; - walk crhmc_walk = walk(problem, p, input.df, input.f, params); + walk crhmc_walk = walk(problem, p, input.df, input.f, params); typedef CrhmcRandomPointGenerator RandomPointGenerator; - RandomPointGenerator::apply(problem, p, nburns, walk_len, randPoints, + RandomPointGenerator::apply(problem, p, nburns, walk_len, randPoints, push_back_policy, rng, F, f, params, crhmc_walk); //crhmc_walk.disable_adaptive(); randPoints.clear(); - RandomPointGenerator::apply(problem, p, rnum, walk_len, randPoints, + RandomPointGenerator::apply(problem, p, rnum, walk_len, randPoints, push_back_policy, rng, F, f, params, crhmc_walk, simdLen, raw_output); } #include "ode_solvers/ode_solvers.hpp" @@ -437,7 +437,7 @@ template < typename CRHMCWalk, int simdLen=1 > -void execute_crhmc(Polytope &P, RNGType &rng, PointList &randPoints, +void execute_crhmc(Polytope &P, RNGType &rng, PointList &randPoints, unsigned int const& walkL, unsigned int const& numpoints, unsigned int const& nburns, NegativeGradientFunctor *F=NULL, NegativeLogprobFunctor *f=NULL, HessianFunctor *h=NULL, bool raw_output= false){ @@ -471,7 +471,7 @@ crhmc_sampling < NegativeGradientFunctor, simdLen > ->(randPoints, P, rng, walkL, numpoints, nburns, *F, *f, *h, simdLen, raw_output); +>(randPoints, P, rng, walkL, numpoints, nburns, *F, *f, *h, simdLen, raw_output); }else{ typedef crhmc_input < @@ -500,7 +500,7 @@ crhmc_sampling < NegativeGradientFunctor, simdLen > ->(randPoints, P, rng, walkL, numpoints, nburns, *F, *f, zerof, simdLen, raw_output); +>(randPoints, P, rng, walkL, numpoints, nburns, *F, *f, zerof, simdLen, raw_output); } } template