From 7c967086f2bbbc0b4c62fbd82c54ef915b86fd9c Mon Sep 17 00:00:00 2001 From: akisschinas Date: Tue, 15 Oct 2024 11:07:38 +0300 Subject: [PATCH] Move mmcs functionality to include from examples --- examples/mmcs_method/CMakeLists.txt | 12 +- ...hpoly_50_dim_parallel.cpp => parallel.cpp} | 0 examples/mmcs_method/random_hpoly_50_dim.cpp | 170 ------------------ examples/mmcs_method/simple.cpp | 63 +++++++ examples/mmcs_method/skinny_cube_10_dim.cpp | 169 ----------------- include/sampling/mmcs.hpp | 149 +++++++++++++++ 6 files changed, 217 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..a614be2a7 100644 --- a/examples/mmcs_method/CMakeLists.txt +++ b/examples/mmcs_method/CMakeLists.txt @@ -108,13 +108,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..d45af4a8e --- /dev/null +++ b/examples/mmcs_method/simple.cpp @@ -0,0 +1,63 @@ +// 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" + + +template +void run_main() +{ + typedef Cartesian Kernel; + typedef boost::mt19937 PolyRNGType; + typedef typename Kernel::Point Point; + typedef HPolytope Hpolytope; + typedef Eigen::Matrix VT; + typedef Eigen::Matrix MT; + + { + int n = 50; + Hpolytope P = random_hpoly(n, 4*n, 127); // we fix the example polytope, seed = 127 + + MT S; + int total_neff; + mmcs(P, 400, S, total_neff); + + 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 n = 10; + Hpolytope P = generate_skinny_cube(n, false); + + MT S; + int total_neff; + mmcs(P, 1000, S, total_neff); + + 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/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..830305b22 100644 --- a/include/sampling/mmcs.hpp +++ b/include/sampling/mmcs.hpp @@ -125,4 +125,153 @@ bool perform_mmcs_step(Polytope &P, return false; } + +template +< + typename Polytope, + typename MT +> +void mmcs(Polytope& P, int const& Neff, MT& S, int& total_neff) +{ + using NT = typename Polytope::NT; + using VT = typename Polytope::VT; + using Point = typename Polytope::PointType; + + const int n = P.dimension(); + MT T = MT::Identity(n, n); + VT T_shift = VT::Zero(n); + + using RNGType = BoostRandomNumberGenerator; + RNGType rng(n); + + unsigned int current_Neff = Neff; + unsigned int round_it = 1; + unsigned int num_rounding_steps = 20 * n; + unsigned int walk_length = 1; + 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