diff --git a/dingo/bindings/bindings.cpp b/dingo/bindings/bindings.cpp index 0e3a62f..2131247 100644 --- a/dingo/bindings/bindings.cpp +++ b/dingo/bindings/bindings.cpp @@ -14,7 +14,7 @@ #include #include "bindings.h" #include "hmc_sampling.h" - +#include "sampling/mmcs.hpp" using namespace std; @@ -134,13 +134,19 @@ double HPolytopeCPP::apply_sampling(int walk_len, } else if (strcmp(method, "vaidya_walk")) { // vaidya walk uniform_sampling(rand_points, HP, rng, walk_len, number_of_points, starting_point, number_of_points_to_burn); - } else if (strcmp(method, "mmcs")) { // vaidya walk - MT S; - int total_ess; - //TODO: avoid passing polytopes as non-const references - const Hpolytope HP_const = HP; - mmcs(HP_const, ess, S, total_ess, walk_len, rng); - samples = S.data(); + } else if (strcmp(method, "mmcs") == 0) { + // Use volesti's MMCS implementation + MT S; + int total_neff = 0; + mmcs(HP, ess, S, total_neff, walk_len, rng); + + // Copy results to output array + for (int i = 0; i < d; i++) { + for (int j = 0; j < S.cols(); j++) { + samples[i * S.cols() + j] = S(i, j); + } + } + return total_neff; } else if (strcmp(method, "gaussian_hmc_walk")) { // Gaussian sampling with exact HMC walk NT a = NT(1)/(NT(2)*variance); gaussian_sampling(rand_points, HP, rng, walk_len, number_of_points, a, @@ -202,313 +208,5 @@ void HPolytopeCPP::get_polytope_as_matrices(double* new_A, double* new_b) const for (int i=0; i < n_hyperplanes; i++){ new_b[i] = new_b_temp[i]; } - -} - - -void HPolytopeCPP::mmcs_initialize(int d, int ess, bool psrf_check, bool parallelism, int num_threads) { - - mmcs_set_of_parameters = mmcs_params(d, ess, psrf_check, parallelism, num_threads); - } - -double HPolytopeCPP::mmcs_step(double* inner_point, double radius, int &N) { - - HP.normalize(); - int d = HP.dimension(); - - VT inner_vec(d); - NT max_s; - - for (int i = 0; i < d; i++){ - inner_vec(i) = inner_point[i]; - } - - Point inner_point2(inner_vec); - CheBall = std::pair(inner_point2, radius); - - HP.set_InnerBall(CheBall); - - RNGType rng(d); - - - if (mmcs_set_of_parameters.request_rounding && mmcs_set_of_parameters.rounding_completed) - { - mmcs_set_of_parameters.req_round_temp = false; - } - - if (mmcs_set_of_parameters.req_round_temp) - { - mmcs_set_of_parameters.nburns = mmcs_set_of_parameters.num_rounding_steps / mmcs_set_of_parameters.window + 1; - } - else - { - mmcs_set_of_parameters.nburns = mmcs_set_of_parameters.max_num_samples / mmcs_set_of_parameters.window + 1; - } - - NT L = NT(6) * std::sqrt(NT(d)) * CheBall.second; - AcceleratedBilliardWalk WalkType(L); - - unsigned int Neff_sampled; - MT TotalRandPoints; - if (mmcs_set_of_parameters.parallelism) - { - mmcs_set_of_parameters.complete = perform_parallel_mmcs_step(HP, rng, mmcs_set_of_parameters.walk_length, - mmcs_set_of_parameters.Neff, - mmcs_set_of_parameters.max_num_samples, - mmcs_set_of_parameters.window, - Neff_sampled, mmcs_set_of_parameters.total_samples, - mmcs_set_of_parameters.num_rounding_steps, - TotalRandPoints, CheBall.first, mmcs_set_of_parameters.nburns, - mmcs_set_of_parameters.num_threads, - mmcs_set_of_parameters.req_round_temp, L); - } - else - { - mmcs_set_of_parameters.complete = perform_mmcs_step(HP, rng, mmcs_set_of_parameters.walk_length, mmcs_set_of_parameters.Neff, - mmcs_set_of_parameters.max_num_samples, mmcs_set_of_parameters.window, - Neff_sampled, mmcs_set_of_parameters.total_samples, - mmcs_set_of_parameters.num_rounding_steps, TotalRandPoints, CheBall.first, - mmcs_set_of_parameters.nburns, mmcs_set_of_parameters.req_round_temp, WalkType); - } - - mmcs_set_of_parameters.store_ess(mmcs_set_of_parameters.phase) = Neff_sampled; - mmcs_set_of_parameters.store_nsamples(mmcs_set_of_parameters.phase) = mmcs_set_of_parameters.total_samples; - mmcs_set_of_parameters.phase++; - mmcs_set_of_parameters.Neff -= Neff_sampled; - std::cout << "phase " << mmcs_set_of_parameters.phase << ": number of correlated samples = " << mmcs_set_of_parameters.total_samples << ", effective sample size = " << Neff_sampled; - mmcs_set_of_parameters.total_neff += Neff_sampled; - - mmcs_set_of_parameters.samples.conservativeResize(d, mmcs_set_of_parameters.total_number_of_samples_in_P0 + mmcs_set_of_parameters.total_samples); - for (int i = 0; i < mmcs_set_of_parameters.total_samples; i++) - { - mmcs_set_of_parameters.samples.col(i + mmcs_set_of_parameters.total_number_of_samples_in_P0) = - mmcs_set_of_parameters.T * TotalRandPoints.row(i).transpose() + mmcs_set_of_parameters.T_shift; - } - - N = mmcs_set_of_parameters.total_number_of_samples_in_P0 + mmcs_set_of_parameters.total_samples; - mmcs_set_of_parameters.total_number_of_samples_in_P0 += mmcs_set_of_parameters.total_samples; - - if (!mmcs_set_of_parameters.complete) - { - if (mmcs_set_of_parameters.request_rounding && !mmcs_set_of_parameters.rounding_completed) - { - VT shift(d), s(d); - MT V(d, d), S(d, d), round_mat; - for (int i = 0; i < d; ++i) - { - shift(i) = TotalRandPoints.col(i).mean(); - } - - for (int i = 0; i < mmcs_set_of_parameters.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(d); - V = MT::Identity(d, d); - } - max_s = s.maxCoeff(); - S = s.asDiagonal(); - round_mat = V * S; - - mmcs_set_of_parameters.round_it++; - HP.shift(shift); - HP.linear_transformIt(round_mat); - mmcs_set_of_parameters.T_shift += mmcs_set_of_parameters.T * shift; - mmcs_set_of_parameters.T = mmcs_set_of_parameters.T * round_mat; - - std::cout << ", ratio of the maximum singilar value over the minimum singular value = " << max_s << std::endl; - - if (max_s <= mmcs_set_of_parameters.s_cutoff || mmcs_set_of_parameters.round_it > mmcs_set_of_parameters.num_its) - { - mmcs_set_of_parameters.rounding_completed = true; - } - } - else - { - std::cout<<"\n"; - } - } - else if (!mmcs_set_of_parameters.psrf_check) - { - NT max_psrf = univariate_psrf(mmcs_set_of_parameters.samples).maxCoeff(); - std::cout << "\n[5]total ess " << mmcs_set_of_parameters.total_neff << ": number of correlated samples = " << mmcs_set_of_parameters.samples.cols()<(mmcs_set_of_parameters.samples).maxCoeff() << std::endl; - std::cout<<"\n\n"; - return 1.5; - } - else - { - TotalRandPoints.resize(0, 0); - NT max_psrf = univariate_psrf(mmcs_set_of_parameters.samples).maxCoeff(); - - if (max_psrf < 1.1 && mmcs_set_of_parameters.total_neff >= mmcs_set_of_parameters.fixed_Neff) { - std::cout << "\n[4]total ess " << mmcs_set_of_parameters.total_neff << ": number of correlated samples = " << mmcs_set_of_parameters.samples.cols()<(mmcs_set_of_parameters.samples).maxCoeff() << std::endl; - std::cout<<"\n\n"; - return 1.5; - } - std::cerr << "\n [1]maximum marginal PSRF: " << max_psrf << std::endl; - - while (max_psrf > 1.1 && mmcs_set_of_parameters.total_neff >= mmcs_set_of_parameters.fixed_Neff) { - - mmcs_set_of_parameters.Neff += mmcs_set_of_parameters.store_ess(mmcs_set_of_parameters.skip_phase); - mmcs_set_of_parameters.total_neff -= mmcs_set_of_parameters.store_ess(mmcs_set_of_parameters.skip_phase); - - mmcs_set_of_parameters.total_number_of_samples_in_P0 -= mmcs_set_of_parameters.store_nsamples(mmcs_set_of_parameters.skip_phase); - N -= mmcs_set_of_parameters.store_nsamples(mmcs_set_of_parameters.skip_phase); - - MT S = mmcs_set_of_parameters.samples; - mmcs_set_of_parameters.samples.resize(d, mmcs_set_of_parameters.total_number_of_samples_in_P0); - mmcs_set_of_parameters.samples = - S.block(0, mmcs_set_of_parameters.store_nsamples(mmcs_set_of_parameters.skip_phase), d, mmcs_set_of_parameters.total_number_of_samples_in_P0); - - mmcs_set_of_parameters.skip_phase++; - - max_psrf = univariate_psrf(mmcs_set_of_parameters.samples).maxCoeff(); - - std::cerr << "[2]maximum marginal PSRF: " << max_psrf << std::endl; - std::cerr << "[2]total ess: " << mmcs_set_of_parameters.total_neff << std::endl; - - if (max_psrf < 1.1 && mmcs_set_of_parameters.total_neff >= mmcs_set_of_parameters.fixed_Neff) { - return 1.5; - } - } - std::cout << "[3]total ess " << mmcs_set_of_parameters.total_neff << ": number of correlated samples = " << mmcs_set_of_parameters.samples.cols()<(mmcs_set_of_parameters.samples).maxCoeff() << std::endl; - std::cout<<"\n\n"; - return 0.0; - } - - return 0.0; -} - -void HPolytopeCPP::get_mmcs_samples(double* T_matrix, double* T_shift, double* samples) { - - int n_variables = HP.dimension(); - - int t_mat_index = 0; - for (int i = 0; i < n_variables; i++){ - for (int j = 0; j < n_variables; j++){ - T_matrix[t_mat_index++] = mmcs_set_of_parameters.T(i, j); - } - } - - // create the shift vector - for (int i = 0; i < n_variables; i++){ - T_shift[i] = mmcs_set_of_parameters.T_shift[i]; - } - - int N = mmcs_set_of_parameters.samples.cols(); - - int t_si = 0; - for (int i = 0; i < n_variables; i++){ - for (int j = 0; j < N; j++){ - samples[t_si++] = mmcs_set_of_parameters.samples(i, j); - } - } - mmcs_set_of_parameters.samples.resize(0,0); -} - - -////////// Start of "rounding()" ////////// -void HPolytopeCPP::apply_rounding(int rounding_method, double* new_A, double* new_b, - double* T_matrix, double* shift, double &round_value, - double* inner_point, double radius){ - - // make a copy of the initial HP which will be used for the rounding step - auto P(HP); - RNGType rng(P.dimension()); - P.normalize(); - - - - // read the inner point provided by the user and the radius - int d = P.dimension(); - VT inner_vec(d); - - for (int i = 0; i < d; i++){ - inner_vec(i) = inner_point[i]; - } - - Point inner_point2(inner_vec); - CheBall = std::pair(inner_point2, radius); - P.set_InnerBall(CheBall); - - // set the output variable of the rounding step - round_result round_res; - - // walk length will always be equal to 2 - int walk_len = 2; - - // run the rounding method - if (rounding_method == 1) { // max ellipsoid - round_res = inscribed_ellipsoid_rounding(P, CheBall.first); - - } else if (rounding_method == 2) { // isotropization - round_res = svd_rounding(P, CheBall, 1, rng); - } else if (rounding_method == 3) { // min ellipsoid - round_res = min_sampling_covering_ellipsoid_rounding(P, - CheBall, - walk_len, - rng); - } else { - throw std::runtime_error("Unknown rounding method."); - } - - // create the new_A matrix - MT A_to_copy = P.get_mat(); - int n_hyperplanes = P.num_of_hyperplanes(); - int n_variables = P.dimension(); - - auto n_si = 0; - for (int i = 0; i < n_hyperplanes; i++){ - for (int j = 0; j < n_variables; j++){ - new_A[n_si++] = A_to_copy(i,j); - } - } - - // create the new_b vector - VT new_b_temp = P.get_vec(); - for (int i=0; i < n_hyperplanes; i++){ - new_b[i] = new_b_temp[i]; - } - - // create the T matrix - MT T_matrix_temp = get<0>(round_res); - auto t_si = 0; - for (int i = 0; i < n_variables; i++){ - for (int j = 0; j < n_variables; j++){ - T_matrix[t_si++] = T_matrix_temp(i,j); - } - } - - // create the shift vector - VT shift_temp = get<1>(round_res); - for (int i = 0; i < n_variables; i++){ - shift[i] = shift_temp[i]; - } - - // get the round val value from the rounding step and print int - round_value = get<2>(round_res); - -} -////////// End of "rounding()" ////////// diff --git a/dingo/bindings/bindings.h b/dingo/bindings/bindings.h index 553ea99..18c8af0 100644 --- a/dingo/bindings/bindings.h +++ b/dingo/bindings/bindings.h @@ -48,72 +48,6 @@ typedef typename Hpolytope::MT MT; typedef typename Hpolytope::VT VT; typedef BoostRandomNumberGenerator RNGType; - -template -struct mmcs_parameters -{ -public: - - mmcs_parameters() {} - - mmcs_parameters(int d, int ess, bool _psrf_check, bool _parallelism, int _num_threads) - : T(MT::Identity(d,d)) - , T_shift(VT::Zero(d)) - , store_ess(VT::Zero(50)) - , store_nsamples(VT::Zero(50)) - , skip_phase(0) - , num_rounding_steps(20*d) - , walk_length(1) - , num_its(20) - , Neff(ess) - , fixed_Neff(ess) - , phase(0) - , window(100) - , max_num_samples(100 * d) - , round_it(1) - , total_number_of_samples_in_P0(0) - , total_neff(0) - , num_threads(_num_threads) - , psrf_check(_psrf_check) - , parallelism(_parallelism) - , complete(false) - , request_rounding(true) - , rounding_completed(false) - , s_cutoff(NT(3)) - { - req_round_temp = request_rounding; - } - - MT T; - MT samples; - VT T_shift; - VT store_ess; - VT store_nsamples; - unsigned int skip_phase; - unsigned int num_rounding_steps; - unsigned int walk_length; - unsigned int num_its; - int Neff; - int fixed_Neff; - unsigned int phase; - unsigned int window; - unsigned int max_num_samples; - unsigned int total_samples; - unsigned int nburns; - unsigned int round_it; - unsigned int total_number_of_samples_in_P0; - unsigned int total_neff; - unsigned int num_threads; - bool psrf_check; - bool parallelism; - bool complete; - bool request_rounding; - bool rounding_completed; - bool req_round_temp; - NT s_cutoff; -}; - - // This is the HPolytopeCPP class; the main volesti class that is running the compute_volume(), rounding() and sampling() methods class HPolytopeCPP{ @@ -121,12 +55,6 @@ class HPolytopeCPP{ std::pair CheBall; - // regarding the rounding step - typedef std::tuple round_result; - typedef mmcs_parameters mmcs_params; - - mmcs_params mmcs_set_of_parameters; - // The class and its main specs HPolytopeCPP(); HPolytopeCPP(double *A, double *b, int n_hyperplanes, int n_variables); @@ -143,19 +71,7 @@ class HPolytopeCPP{ char* method, double* inner_point, double radius, double* samples, double variance_value, double* bias_vector, int ess); - void mmcs_initialize(int d, int ess, bool psrf_check, bool parallelism, int num_threads); - - double mmcs_step(double* inner_point_for_c, double radius, int &N); - - void get_mmcs_samples(double* T_matrix, double* T_shift, double* samples); - void get_polytope_as_matrices(double* new_A, double* new_b) const; - - // the rounding() function - void apply_rounding(int rounding_method, double* new_A, double* new_b, double* T_matrix, - double* shift, double &round_value, double* inner_point, double radius); - }; - #endif