From 2626ef4874bdac20103568c8d71fbeb521c680b8 Mon Sep 17 00:00:00 2001 From: akisschinas Date: Tue, 15 Oct 2024 11:07:38 +0300 Subject: [PATCH] Create a mmcs sampling function --- examples/mmcs_method/CMakeLists.txt | 13 +- ...hpoly_50_dim_parallel.cpp => parallel.cpp} | 0 examples/mmcs_method/random_hpoly_50_dim.cpp | 170 ----------------- examples/mmcs_method/simple.cpp | 67 +++++++ examples/mmcs_method/skinny_cube_10_dim.cpp | 169 ---------------- include/sampling/mmcs.hpp | 180 ++++++++++++++++++ 6 files changed, 253 insertions(+), 346 deletions(-) rename examples/mmcs_method/{random_hpoly_50_dim_parallel.cpp => parallel.cpp} (100%) delete mode 100644 examples/mmcs_method/random_hpoly_50_dim.cpp create mode 100644 examples/mmcs_method/simple.cpp delete mode 100644 examples/mmcs_method/skinny_cube_10_dim.cpp diff --git a/examples/mmcs_method/CMakeLists.txt b/examples/mmcs_method/CMakeLists.txt index 5ecd129f2..06e0137ca 100644 --- a/examples/mmcs_method/CMakeLists.txt +++ b/examples/mmcs_method/CMakeLists.txt @@ -99,6 +99,7 @@ else () add_definitions(${CMAKE_CXX_FLAGS} "-std=c++17") # enable C++17 standard add_definitions(${CMAKE_CXX_FLAGS} "-O3") # optimization of the compiler + add_definitions(${CMAKE_CXX_FLAGS} "-g") #add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lgsl") add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lm") add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-ldl") @@ -108,13 +109,11 @@ else () find_package(OpenMP REQUIRED) - add_executable (skinny_cube_10_dim skinny_cube_10_dim.cpp) - add_executable (random_hpoly_50_dim random_hpoly_50_dim.cpp) - add_executable (random_hpoly_50_dim_parallel random_hpoly_50_dim_parallel.cpp) - - TARGET_LINK_LIBRARIES(skinny_cube_10_dim ${LP_SOLVE}) - TARGET_LINK_LIBRARIES(random_hpoly_50_dim ${LP_SOLVE}) - TARGET_LINK_LIBRARIES(random_hpoly_50_dim_parallel ${LP_SOLVE} OpenMP::OpenMP_CXX) + add_executable (simple simple.cpp) + add_executable (parallel parallel.cpp) + TARGET_LINK_LIBRARIES(simple ${LP_SOLVE}) + TARGET_LINK_LIBRARIES(parallel ${LP_SOLVE}) + TARGET_LINK_LIBRARIES(parallel ${LP_SOLVE} OpenMP::OpenMP_CXX) endif() diff --git a/examples/mmcs_method/random_hpoly_50_dim_parallel.cpp b/examples/mmcs_method/parallel.cpp similarity index 100% rename from examples/mmcs_method/random_hpoly_50_dim_parallel.cpp rename to examples/mmcs_method/parallel.cpp diff --git a/examples/mmcs_method/random_hpoly_50_dim.cpp b/examples/mmcs_method/random_hpoly_50_dim.cpp deleted file mode 100644 index 90cae4fa3..000000000 --- a/examples/mmcs_method/random_hpoly_50_dim.cpp +++ /dev/null @@ -1,170 +0,0 @@ -// VolEsti (volume computation and sampling library) - -// Copyright (c) 2012-2021 Vissarion Fisikopoulos -// Copyright (c) 2018-2021 Apostolos Chalkis - -#include "Eigen/Eigen" - -#include -#include -#include -#include -#include -#include "random_walks/random_walks.hpp" -#include "volume/volume_sequence_of_balls.hpp" -#include "volume/volume_cooling_gaussians.hpp" -#include "sampling/mmcs.hpp" -#include "generators/h_polytopes_generator.h" -#include "diagnostics/multivariate_psrf.hpp" -#include "diagnostics/univariate_psrf.hpp" -#include "diagnostics/ess_window_updater.hpp" - - -template -void run_main() -{ - typedef Cartesian Kernel; - typedef BoostRandomNumberGenerator RNGType; - typedef boost::mt19937 PolyRNGType; - typedef typename Kernel::Point Point; - typedef HPolytope Hpolytope; - typedef Eigen::Matrix VT; - typedef Eigen::Matrix MT; - - int n = 50; - - RNGType rng(n); - Hpolytope P = random_hpoly(n, 4*n, 127); // we fix the example polytope, seed = 127 - std::list randPoints; - - MT T = MT::Identity(n, n); - VT T_shift = VT::Zero(n); - - unsigned int round_it = 1, num_rounding_steps = 20*n, - walk_length = 1, num_its = 20, Neff = 4000, total_neff = 0, phase = 0, - window = 100, max_num_samples = 100 * n, total_samples, nburns, total_number_of_samples_in_P0 = 0; - NT max_s, s_cutoff = 3.0, L; - bool complete = false, request_rounding = true, - rounding_completed = false; - bool req_round_temp = request_rounding; - - std::pair InnerBall; - - MT S; - - std::cout << "target effective sample size = " << Neff << "\n" << std::endl; - - while(true) - { - phase++; - if (request_rounding && rounding_completed) - { - req_round_temp = false; - } - - if (req_round_temp) - { - nburns = num_rounding_steps / window + 1; - } - else - { - nburns = max_num_samples / window + 1; - } - - InnerBall = P.ComputeInnerBall(); - L = NT(6) * std::sqrt(NT(n)) * InnerBall.second; - AcceleratedBilliardWalk WalkType(L); - - unsigned int Neff_sampled; - MT TotalRandPoints; - complete = perform_mmcs_step(P, rng, walk_length, Neff, max_num_samples, window, - Neff_sampled, total_samples, num_rounding_steps, TotalRandPoints, - InnerBall.first, nburns, req_round_temp, WalkType); - - Neff -= Neff_sampled; - std::cout << "phase " << phase << ": number of correlated samples = " << total_samples << ", effective sample size = " << Neff_sampled; - total_neff += Neff_sampled; - Neff_sampled = 0; - - MT Samples = TotalRandPoints.transpose(); //do not copy TODO! - for (int i = 0; i < total_samples; i++) - { - Samples.col(i) = T * Samples.col(i) + T_shift; - } - - S.conservativeResize(P.dimension(), total_number_of_samples_in_P0 + total_samples); - S.block(0, total_number_of_samples_in_P0, P.dimension(), total_samples) = Samples.block(0, 0, P.dimension(), total_samples); - total_number_of_samples_in_P0 += total_samples; - if (!complete) - { - if (request_rounding && !rounding_completed) - { - VT shift(n), s(n); - MT V(n,n), S(n,n), round_mat; - for (int i = 0; i < P.dimension(); ++i) - { - shift(i) = TotalRandPoints.col(i).mean(); - } - - for (int i = 0; i < total_samples; ++i) - { - TotalRandPoints.row(i) = TotalRandPoints.row(i) - shift.transpose(); - } - - Eigen::BDCSVD svd(TotalRandPoints, Eigen::ComputeFullV); - s = svd.singularValues() / svd.singularValues().minCoeff(); - - if (s.maxCoeff() >= 2.0) - { - for (int i = 0; i < s.size(); ++i) - { - if (s(i) < 2.0) - { - s(i) = 1.0; - } - } - V = svd.matrixV(); - } - else - { - s = VT::Ones(P.dimension()); - V = MT::Identity(P.dimension(), P.dimension()); - } - max_s = s.maxCoeff(); - S = s.asDiagonal(); - round_mat = V * S; - - round_it++; - P.shift(shift); - P.linear_transformIt(round_mat); - T_shift += T * shift; - T = T * round_mat; - - std::cout << ", ratio of the maximum singilar value over the minimum singular value = " << max_s << std::endl; - - if (max_s <= s_cutoff || round_it > num_its) - { - rounding_completed = true; - } - } - else - { - std::cout<<"\n"; - } - } - else - { - std::cout<<"\n\n"; - break; - } - } - - std::cerr << "sum of effective sample sizes: " << total_neff << std::endl; - std::cerr << "multivariate PSRF: " << multivariate_psrf(S) << std::endl; - std::cerr << "maximum marginal PSRF: " << univariate_psrf(S).maxCoeff() << std::endl; -} - -int main() { - run_main(); - return 0; -} diff --git a/examples/mmcs_method/simple.cpp b/examples/mmcs_method/simple.cpp new file mode 100644 index 000000000..cd11d4c69 --- /dev/null +++ b/examples/mmcs_method/simple.cpp @@ -0,0 +1,67 @@ +// VolEsti (volume computation and sampling library) + +// Copyright (c) 2012-2021 Vissarion Fisikopoulos +// Copyright (c) 2018-2021 Apostolos Chalkis + +#include "Eigen/Eigen" + +#include +#include +#include +#include +#include +#include "random_walks/random_walks.hpp" +#include "volume/volume_sequence_of_balls.hpp" +#include "volume/volume_cooling_gaussians.hpp" +#include "sampling/mmcs.hpp" +#include "generators/h_polytopes_generator.h" +#include "generators/known_polytope_generators.h" +#include "diagnostics/multivariate_psrf.hpp" +#include "diagnostics/univariate_psrf.hpp" +#include "diagnostics/ess_window_updater.hpp" + +typedef double NT; +typedef Cartesian Kernel; +typedef typename Kernel::Point Point; +typedef HPolytope Hpolytope; + +template +void mmcs_sampling(Polytope const& P) +{ + typedef Eigen::Matrix VT; + typedef Eigen::Matrix MT; + + MT S; + int total_neff; + + mmcs(P, 100, S, total_neff); + + for (int i = 0; i < S.cols(); ++i) { + Point p(S.rows()); + for (int j = 0; j < S.rows(); ++j) { + p.set_coord(j, S(j,i)); + } + if (P.is_in(p) == 0) { + std::cout << "Sample point out of the polytope."; + } + } + + std::cerr << "sum of effective sample sizes: " << total_neff << std::endl; + std::cerr << "multivariate PSRF: " << multivariate_psrf(S) << std::endl; + std::cerr << "maximum marginal PSRF: " << univariate_psrf(S).maxCoeff() << std::endl; +} + +int main() +{ + typedef boost::mt19937 PolyRNGType; + + int n = 50; + Hpolytope P1 = random_hpoly(n, 4*n, 127); // we fix the example polytope, seed = 127 + mmcs_sampling(P1); + + Hpolytope P2 = generate_skinny_cube(10, false); + mmcs_sampling(P2); + + return 0; +} + diff --git a/examples/mmcs_method/skinny_cube_10_dim.cpp b/examples/mmcs_method/skinny_cube_10_dim.cpp deleted file mode 100644 index 036b10a34..000000000 --- a/examples/mmcs_method/skinny_cube_10_dim.cpp +++ /dev/null @@ -1,169 +0,0 @@ -// VolEsti (volume computation and sampling library) - -// Copyright (c) 2012-2021 Vissarion Fisikopoulos -// Copyright (c) 2018-2021 Apostolos Chalkis - -#include "Eigen/Eigen" - -#include -#include -#include -#include -#include -#include "random_walks/random_walks.hpp" -#include "volume/volume_sequence_of_balls.hpp" -#include "volume/volume_cooling_gaussians.hpp" -#include "sampling/mmcs.hpp" -#include "generators/known_polytope_generators.h" -#include "diagnostics/multivariate_psrf.hpp" -#include "diagnostics/univariate_psrf.hpp" -#include "diagnostics/ess_window_updater.hpp" - - -template -void run_main() -{ - typedef Cartesian Kernel; - typedef BoostRandomNumberGenerator RNGType; - typedef typename Kernel::Point Point; - typedef HPolytope Hpolytope; - typedef Eigen::Matrix VT; - typedef Eigen::Matrix MT; - - int n = 10; - - Hpolytope P = generate_skinny_cube(n, false); - std::list randPoints; - RNGType rng(n); - - MT T = MT::Identity(n, n); - VT T_shift = VT::Zero(n); - - unsigned int round_it = 1, num_rounding_steps = 20*n, - walk_length = 1, num_its = 20, Neff = 1000, total_neff = 0, phase = 0, - window = 100, max_num_samples = 100 * n, total_samples, nburns, total_number_of_samples_in_P0 = 0; - NT max_s, s_cutoff = 3.0, L; - bool complete = false, request_rounding = true, - rounding_completed = false; - bool req_round_temp = request_rounding; - - std::pair InnerBall; - - MT S; - - std::cout << "target effective sample size = " << Neff << "\n" << std::endl; - - while(true) - { - phase++; - if (request_rounding && rounding_completed) - { - req_round_temp = false; - } - - if (req_round_temp) - { - nburns = num_rounding_steps / window + 1; - } - else - { - nburns = max_num_samples / window + 1; - } - - InnerBall = P.ComputeInnerBall(); - L = NT(6) * std::sqrt(NT(n)) * InnerBall.second; - AcceleratedBilliardWalk WalkType(L); - - unsigned int Neff_sampled; - MT TotalRandPoints; - complete = perform_mmcs_step(P, rng, walk_length, Neff, max_num_samples, window, - Neff_sampled, total_samples, num_rounding_steps, TotalRandPoints, - InnerBall.first, nburns, req_round_temp, WalkType); - - Neff -= Neff_sampled; - std::cout << "phase " << phase << ": number of correlated samples = " << total_samples << ", effective sample size = " << Neff_sampled; - total_neff += Neff_sampled; - Neff_sampled = 0; - - MT Samples = TotalRandPoints.transpose(); //do not copy TODO! - for (int i = 0; i < total_samples; i++) - { - Samples.col(i) = T * Samples.col(i) + T_shift; - } - - S.conservativeResize(P.dimension(), total_number_of_samples_in_P0 + total_samples); - S.block(0, total_number_of_samples_in_P0, P.dimension(), total_samples) = Samples.block(0, 0, P.dimension(), total_samples); - total_number_of_samples_in_P0 += total_samples; - if (!complete) - { - if (request_rounding && !rounding_completed) - { - VT shift(n), s(n); - MT V(n,n), S(n,n), round_mat; - for (int i = 0; i < P.dimension(); ++i) - { - shift(i) = TotalRandPoints.col(i).mean(); - } - - for (int i = 0; i < total_samples; ++i) - { - TotalRandPoints.row(i) = TotalRandPoints.row(i) - shift.transpose(); - } - - Eigen::BDCSVD svd(TotalRandPoints, Eigen::ComputeFullV); - s = svd.singularValues() / svd.singularValues().minCoeff(); - - if (s.maxCoeff() >= 2.0) - { - for (int i = 0; i < s.size(); ++i) - { - if (s(i) < 2.0) - { - s(i) = 1.0; - } - } - V = svd.matrixV(); - } - else - { - s = VT::Ones(P.dimension()); - V = MT::Identity(P.dimension(), P.dimension()); - } - max_s = s.maxCoeff(); - S = s.asDiagonal(); - round_mat = V * S; - - round_it++; - P.shift(shift); - P.linear_transformIt(round_mat); - T_shift += T * shift; - T = T * round_mat; - - std::cout << ", ratio of the maximum singilar value over the minimum singular value = " << max_s << std::endl; - - if (max_s <= s_cutoff || round_it > num_its) - { - rounding_completed = true; - } - } - else - { - std::cout<<"\n"; - } - } - else - { - std::cout<<"\n\n"; - break; - } - } - - std::cerr << "sum of effective sample sizes: " << total_neff << std::endl; - std::cerr << "multivariate PSRF: " << multivariate_psrf(S) << std::endl; - std::cerr << "maximum marginal PSRF: " << univariate_psrf(S).maxCoeff() << std::endl; -} - -int main() { - run_main(); - return 0; -} diff --git a/include/sampling/mmcs.hpp b/include/sampling/mmcs.hpp index 02791535b..0e8cc990f 100644 --- a/include/sampling/mmcs.hpp +++ b/include/sampling/mmcs.hpp @@ -125,4 +125,184 @@ bool perform_mmcs_step(Polytope &P, return false; } +template +< + typename Polytope, + typename MT +> +void mmcs(Polytope const& Pin, + int const& Neff, + MT& S, + int& total_neff) +{ + mmcs(Pin, Neff, S, total_neff, 1); +} + +template +< + typename Polytope, + typename MT +> +void mmcs(Polytope const& Pin, + int const& Neff, + MT& S, + int& total_neff, + unsigned int const& walk_length) +{ + using RNGType = BoostRandomNumberGenerator; + RNGType rng(Pin.dimension()); + mmcs(Pin, Neff, S, total_neff, walk_length, rng); +} + +template +< + typename Polytope, + typename MT, + typename RandomNumberGenerator +> +void mmcs(Polytope const& Pin, + int const& Neff, + MT& S, + int& total_neff, + unsigned int const& walk_length, + RandomNumberGenerator &rng) +{ + using NT = typename Polytope::NT; + using VT = typename Polytope::VT; + using Point = typename Polytope::PointType; + + auto P = Pin; + const int n = P.dimension(); + MT T = MT::Identity(n, n); + VT T_shift = VT::Zero(n); + + unsigned int current_Neff = Neff; + unsigned int round_it = 1; + unsigned int num_rounding_steps = 20 * n; + unsigned int num_its = 20; + unsigned int phase = 0; + unsigned int window = 100; + unsigned int max_num_samples = 100 * n; + unsigned int total_samples; + unsigned int nburns; + unsigned int total_number_of_samples_in_P0 = 0; + + NT max_s; + NT s_cutoff = 3.0; + NT L; + bool complete = false; + bool request_rounding = true; + bool rounding_completed = false; + bool req_round_temp = request_rounding; + + std::pair InnerBall; + + std::cout << "target effective sample size = " << current_Neff << "\n" << std::endl; + + while(true) + { + phase++; + if (request_rounding && rounding_completed) + { + req_round_temp = false; + } + + if (req_round_temp) + { + nburns = num_rounding_steps / window + 1; + } + else + { + nburns = max_num_samples / window + 1; + } + + InnerBall = P.ComputeInnerBall(); + L = NT(6) * std::sqrt(NT(n)) * InnerBall.second; + AcceleratedBilliardWalk WalkType(L); + + unsigned int Neff_sampled; + MT TotalRandPoints; + complete = perform_mmcs_step(P, rng, walk_length, current_Neff, max_num_samples, window, + Neff_sampled, total_samples, num_rounding_steps, TotalRandPoints, + InnerBall.first, nburns, req_round_temp, WalkType); + + current_Neff -= Neff_sampled; + std::cout << "phase " << phase << ": number of correlated samples = " << total_samples << ", effective sample size = " << Neff_sampled; + total_neff += Neff_sampled; + Neff_sampled = 0; + + MT Samples = TotalRandPoints.transpose(); //do not copy TODO! + for (int i = 0; i < total_samples; i++) + { + Samples.col(i) = T * Samples.col(i) + T_shift; + } + + S.conservativeResize(P.dimension(), total_number_of_samples_in_P0 + total_samples); + S.block(0, total_number_of_samples_in_P0, P.dimension(), total_samples) = Samples.block(0, 0, P.dimension(), total_samples); + total_number_of_samples_in_P0 += total_samples; + if (!complete) + { + if (request_rounding && !rounding_completed) + { + VT shift(n), s(n); + MT V(n,n), S(n,n), round_mat; + for (int i = 0; i < P.dimension(); ++i) + { + shift(i) = TotalRandPoints.col(i).mean(); + } + + for (int i = 0; i < total_samples; ++i) + { + TotalRandPoints.row(i) = TotalRandPoints.row(i) - shift.transpose(); + } + + Eigen::BDCSVD svd(TotalRandPoints, Eigen::ComputeFullV); + s = svd.singularValues() / svd.singularValues().minCoeff(); + + if (s.maxCoeff() >= 2.0) + { + for (int i = 0; i < s.size(); ++i) + { + if (s(i) < 2.0) + { + s(i) = 1.0; + } + } + V = svd.matrixV(); + } + else + { + s = VT::Ones(P.dimension()); + V = MT::Identity(P.dimension(), P.dimension()); + } + max_s = s.maxCoeff(); + S = s.asDiagonal(); + round_mat = V * S; + + round_it++; + P.shift(shift); + P.linear_transformIt(round_mat); + T_shift += T * shift; + T = T * round_mat; + + std::cout << ", ratio of the maximum singilar value over the minimum singular value = " << max_s << std::endl; + + if (max_s <= s_cutoff || round_it > num_its) + { + rounding_completed = true; + } + } + else + { + std::cout<<"\n"; + } + } + else + { + std::cout<<"\n\n"; + break; + } + } +} + #endif