From 8c1f2b4ad61079b7d4e51d88a2dfe5d08f6d2391 Mon Sep 17 00:00:00 2001 From: Vaibhav Thakkar Date: Wed, 28 Jul 2021 15:34:19 +0530 Subject: [PATCH 1/7] move math helper functions to a common place --- include/volume/math_helpers.h | 44 +++++++++++++++++++++ include/volume/volume_cooling_balls.hpp | 10 +---- include/volume/volume_cooling_gaussians.hpp | 24 +---------- 3 files changed, 47 insertions(+), 31 deletions(-) create mode 100644 include/volume/math_helpers.h diff --git a/include/volume/math_helpers.h b/include/volume/math_helpers.h new file mode 100644 index 000000000..fe1bf8155 --- /dev/null +++ b/include/volume/math_helpers.h @@ -0,0 +1,44 @@ +// VolEsti (volume computation and sampling library) + +// Copyright (c) 2021 Vaibhav Thakkar + +// Contributed and/or modified by Vaibhav Thakkar, as part of Google Summer of Code 2021 program. + +// Licensed under GNU LGPL.3, see LICENCE file + +#ifndef MATH_HELPERS_H +#define MATH_HELPERS_H + +#include +#include + + +//An implementation of Welford's algorithm for mean and variance. +template +std::pair get_mean_variance(std::vector& vec) +{ + NT mean = 0; + NT M2 = 0; + NT variance = 0; + NT delta; + + unsigned int i=0; + for (auto vecit = vec.begin(); vecit!=vec.end(); vecit++, i++) + { + delta = *vecit - mean; + mean += delta / (i + 1); + M2 += delta * (*vecit - mean); + variance = M2 / (i + 1); + } + return std::pair (mean, variance); +} + + +template +static NT log_gamma_function(NT x) +{ + if (x <= NT(100)) return std::log(tgamma(x)); + return (std::log(x - NT(1)) + log_gamma_function(x - NT(1))); +} + +#endif // MATH_HELPERS_H \ No newline at end of file diff --git a/include/volume/volume_cooling_balls.hpp b/include/volume/volume_cooling_balls.hpp index 429655818..f3f41ea34 100644 --- a/include/volume/volume_cooling_balls.hpp +++ b/include/volume/volume_cooling_balls.hpp @@ -22,6 +22,7 @@ #endif #include "convex_bodies/ballintersectconvex.h" #include "sampling/random_point_generators.hpp" +#include "math_helpers.h" //////////////////////////////////// @@ -685,15 +686,6 @@ NT estimate_ratio_interval(PolyBall1 const& Pb1, return NT(ratio_parameters.count_in) / NT(ratio_parameters.tot_count); } - -template -NT log_gamma_function(NT x) -{ - if (x <= NT(100)) return std::log(tgamma(x)); - return (std::log(x - NT(1)) + log_gamma_function(x - NT(1))); -} - - template < typename WalkTypePolicy, diff --git a/include/volume/volume_cooling_gaussians.hpp b/include/volume/volume_cooling_gaussians.hpp index a09ecb24d..aa1e2dd4c 100644 --- a/include/volume/volume_cooling_gaussians.hpp +++ b/include/volume/volume_cooling_gaussians.hpp @@ -23,6 +23,7 @@ #include "random_walks/gaussian_ball_walk.hpp" #include "random_walks/gaussian_cdhr_walk.hpp" #include "sampling/random_point_generators.hpp" +#include "math_helpers.h" /////////////////// Helpers for random walks @@ -54,27 +55,6 @@ struct update_delta> // Springer-Verlag Berlin Heidelberg and The Mathematical Programming Society 2015 // Ben Cousins, Santosh Vempala -//An implementation of Welford's algorithm for mean and variance. -template -std::pair get_mean_variance(std::vector& vec) -{ - NT mean = 0; - NT M2 = 0; - NT variance = 0; - NT delta; - - unsigned int i=0; - for (auto vecit = vec.begin(); vecit!=vec.end(); vecit++, i++) - { - delta = *vecit - mean; - mean += delta / (i + 1); - M2 += delta * (*vecit - mean); - variance = M2 / (i + 1); - } - return std::pair (mean, variance); -} - - // Compute the first variance a_0 for the starting gaussian template void get_first_gaussian(Polytope const& P, @@ -336,7 +316,7 @@ double volume_cooling_gaussians(Polytope const& Pin, // Consider Chebychev center as an internal point auto InnerBall = P.ComputeInnerBall(); if (InnerBall.second < 0.0) return -1.0; - + Point c = InnerBall.first; NT radius = InnerBall.second; From 8cfb29517a9bdf9b40e63fb7a6778a0110955ac5 Mon Sep 17 00:00:00 2001 From: Vaibhav Thakkar Date: Wed, 28 Jul 2021 16:40:14 +0530 Subject: [PATCH 2/7] add utility functions for reading poset and finding an inner point in the order polytope --- include/convex_bodies/orderpolytope.h | 40 +++++++++++++++++ include/misc/misc.h | 65 +++++++++++++++++++++++++++ include/misc/poset.h | 37 +++++++++++++++ 3 files changed, 142 insertions(+) diff --git a/include/convex_bodies/orderpolytope.h b/include/convex_bodies/orderpolytope.h index 31fc01113..33f5c0d6d 100644 --- a/include/convex_bodies/orderpolytope.h +++ b/include/convex_bodies/orderpolytope.h @@ -93,6 +93,18 @@ class OrderPolytope { } + Eigen::SparseMatrix get_mat() const + { + return _A.sparseView(); + } + + + VT get_vec() const + { + return b; + } + + // print polytope in Ax <= b format void print() { @@ -647,6 +659,34 @@ class OrderPolytope { } + // get a point inside the order polytope, + // NOTE: the current implementation only works for non shifted order polytope + VT inner_point() + { + // get topologically sorted list of indices + std::vector sorted_list = poset.topologically_sorted_list(); + + // vector to hold n linearly spaced values between 0-1 + std::vector lin_space_values(_d); + NT start = 0.05, end = 0.95; + NT h = (end - start)/static_cast(_d-1); + NT val = start; + for(auto x=lin_space_values.begin(); x!=lin_space_values.end(); ++x) { + *x = val; + val += h; + } + + // final result vector + VT res(_d); + for(int i=0; i<_d; ++i) { + unsigned int curr_idx = sorted_list[i]; + res(curr_idx) = lin_space_values[i]; + } + + return res; + } + + void normalize() { if (_normalized == true) diff --git a/include/misc/misc.h b/include/misc/misc.h index 7c3c99222..e4578cc7a 100644 --- a/include/misc/misc.h +++ b/include/misc/misc.h @@ -9,6 +9,8 @@ #include #include +#include +#include "poset.h" //function to print rounding to double coordinates template @@ -203,4 +205,67 @@ std::pair read_inner_ball(std::istream &is) { return std::make_pair(center, radius); } +/* read a poset given in the following format: + - First line contains a single positive integer 'n' - number of elements + - Next `m` lines follow containing a pair 'i j' in each line to signify A_i <= A_j + i.e i_th element is less than or equal to the j_th element +*/ +Poset read_poset_from_file(std::string filename) { + typedef typename Poset::RT RT; + typedef typename Poset::RV RV; + + std::ifstream data_file; + data_file.open(filename); + + // read number of elements + unsigned int n; + data_file >> n; + + // read relations line by line + RT curr_relation; + RV relations; + while(data_file >> curr_relation.first >> curr_relation.second) + relations.push_back(curr_relation); + + return Poset(n, relations); +} + + +// read a poset given as an adjacency matrix +std::pair read_poset_from_file_adj_matrix(std::string filename) { + typedef typename Poset::RV RV; + + std::ifstream in; + in.open(filename); + RV edges; + unsigned int x, n = 0; + + // read a single line + std::string line; + std::getline(in, line); + std::stringstream line_ss(line); + while(line_ss >> x) { + if(x) { + edges.emplace_back(0, n); + } + ++n; + } + + // read rest of the lines + for(unsigned int a = 1; a < n; ++a) { + for(unsigned int b = 0; b < n; ++b) { + if(!(in >> x)) { + std::cerr << "Invalid adjacency matrix"; + return std::pair(false, Poset()); + } + + if(x) { + edges.emplace_back(a, b); + } + } + } + + return std::pair(true, Poset(n, edges)); +} + #endif //MISC_H diff --git a/include/misc/poset.h b/include/misc/poset.h index 0f98e3c6c..08296ff18 100644 --- a/include/misc/poset.h +++ b/include/misc/poset.h @@ -13,6 +13,7 @@ #include #include +#include class Poset { @@ -25,6 +26,8 @@ class Poset { RV order_relations; // pairs of form a <= b public: + Poset() {} + Poset(unsigned int _n, RV& _order_relations) : n(_n), order_relations(verify(_order_relations, n)) {} @@ -89,6 +92,40 @@ class Poset { return true; } + + + std::vector topologically_sorted_list() const + { + std::vector > adjList(n); + std::vector indegree(n, 0); + + for(auto x: order_relations) { + adjList[x.first].push_back(x.second); + indegree[x.second] += 1; + } + + std::queue q; + for(unsigned int i=0; i res; + while(!q.empty()) { + unsigned int curr = q.front(); + res.push_back(curr); + q.pop(); + + for(unsigned int i=0; i Date: Fri, 30 Jul 2021 02:21:41 +0530 Subject: [PATCH 3/7] use read poset from misc file --- .../order-polytope-basics/order_polytope.cpp | 19 ++----------------- 1 file changed, 2 insertions(+), 17 deletions(-) diff --git a/examples/order-polytope-basics/order_polytope.cpp b/examples/order-polytope-basics/order_polytope.cpp index 9ea048b2f..c58738f73 100644 --- a/examples/order-polytope-basics/order_polytope.cpp +++ b/examples/order-polytope-basics/order_polytope.cpp @@ -20,21 +20,6 @@ std::vector linspace(T a, T b, size_t N) { typedef typename Poset::RT RT; typedef typename Poset::RV RV; -Poset read_poset_from_file(std::string filename) { - std::ifstream data_file; - data_file.open(filename); - - unsigned int n; - data_file >> n; - - RT curr_relation; - RV relations; - while(data_file >> curr_relation.first >> curr_relation.second) - relations.push_back(curr_relation); - - return Poset(n, relations); -} - int main(int argc, char const *argv[]) { std::cout << "\nPoset operations: \n"; @@ -91,11 +76,11 @@ int main(int argc, char const *argv[]) { ray.set_coord(0, 1.5); std::cout << "incident ray: "; ray.print(); - + OP.compute_reflection(ray, Point(), 2*OP.dimension()); std::cout << "reflected ray: "; ray.print(); - // --------------------------------------------- + // --------------------------------------------- return 0; } From d004421c449e77c00c3d3bf4be2db387201ae2c5 Mon Sep 17 00:00:00 2001 From: Vaibhav Thakkar Date: Tue, 14 Sep 2021 03:36:47 +0530 Subject: [PATCH 4/7] fix header files inclusion --- include/misc/misc.h | 1 + include/random_walks/uniform_billiard_walk.hpp | 1 + include/volume/volume_cooling_balls.hpp | 2 ++ 3 files changed, 4 insertions(+) diff --git a/include/misc/misc.h b/include/misc/misc.h index e4578cc7a..f0a585af3 100644 --- a/include/misc/misc.h +++ b/include/misc/misc.h @@ -10,6 +10,7 @@ #include #include #include +#include #include "poset.h" //function to print rounding to double coordinates diff --git a/include/random_walks/uniform_billiard_walk.hpp b/include/random_walks/uniform_billiard_walk.hpp index 95b542b7d..ca6eb74b7 100644 --- a/include/random_walks/uniform_billiard_walk.hpp +++ b/include/random_walks/uniform_billiard_walk.hpp @@ -22,6 +22,7 @@ #include "convex_bodies/zonoIntersecthpoly.h" #endif #include "sampling/sphere.hpp" +#include "random_walks/boundary_cdhr_walk.hpp" #include "generators/boost_random_number_generator.hpp" #include "sampling/random_point_generators.hpp" #include "volume/sampling_policies.hpp" diff --git a/include/volume/volume_cooling_balls.hpp b/include/volume/volume_cooling_balls.hpp index f3f41ea34..5958f6f98 100644 --- a/include/volume/volume_cooling_balls.hpp +++ b/include/volume/volume_cooling_balls.hpp @@ -20,6 +20,8 @@ #include "convex_bodies/zpolytope.h" #include "convex_bodies/vpolyintersectvpoly.h" #endif +#include "random_walks/uniform_cdhr_walk.hpp" +#include "convex_bodies/ball.h" #include "convex_bodies/ballintersectconvex.h" #include "sampling/random_point_generators.hpp" #include "math_helpers.h" From d53f6dbbb995a75bfe0b6dc7a9eed17afe708a73 Mon Sep 17 00:00:00 2001 From: Vaibhav Thakkar Date: Tue, 14 Sep 2021 03:55:41 +0530 Subject: [PATCH 5/7] add more reformatting changes --- R-proj/src/copula.cpp | 2 +- .../order-polytope-basics/order_polytope.cpp | 5 ++- include/convex_bodies/ellipsoid.h | 32 ++++--------------- include/misc/misc.h | 10 ++---- 4 files changed, 14 insertions(+), 35 deletions(-) diff --git a/R-proj/src/copula.cpp b/R-proj/src/copula.cpp index d2bb8be3a..8b9fd5868 100644 --- a/R-proj/src/copula.cpp +++ b/R-proj/src/copula.cpp @@ -63,7 +63,7 @@ Rcpp::NumericMatrix copula (Rcpp::Nullable r1, typedef boost::mt19937 RNGType; typedef Eigen::Matrix MT; typedef Eigen::Matrix VT; - typedef Ellipsoid CopEll; + typedef Ellipsoid CopEll; unsigned int num_slices = 100, numpoints = 500000; if (m.isNotNull()) { diff --git a/examples/order-polytope-basics/order_polytope.cpp b/examples/order-polytope-basics/order_polytope.cpp index c58738f73..652aca423 100644 --- a/examples/order-polytope-basics/order_polytope.cpp +++ b/examples/order-polytope-basics/order_polytope.cpp @@ -24,7 +24,10 @@ typedef typename Poset::RV RV; int main(int argc, char const *argv[]) { std::cout << "\nPoset operations: \n"; // ----------- basic poset operations ----------- - Poset poset = read_poset_from_file(std::string(argv[1])); + std::string filename (argv[1]); + std::ifstream data_file; + data_file.open(filename); + Poset poset = read_poset_from_file(data_file); poset.print(); std::cout << "Checking if a sample Point (linearly spaced coordinates) lies inside the poset: "; diff --git a/include/convex_bodies/ellipsoid.h b/include/convex_bodies/ellipsoid.h index 3c1b35757..155e6db1f 100644 --- a/include/convex_bodies/ellipsoid.h +++ b/include/convex_bodies/ellipsoid.h @@ -16,23 +16,17 @@ #include #include -#include +#include "volume/math_helpers.h" -template -NT log_gamma_function(NT x) -{ - if (x <= NT(100)) return std::log(tgamma(x)); - return (std::log(x - NT(1)) + log_gamma_function(x - NT(1))); -} - - -template +template class Ellipsoid{ -private: +public: +typedef Point PointType; typedef typename Point::FT NT; typedef typename std::vector::iterator viterator; typedef Eigen::Matrix VT; + typedef Eigen::Matrix MT; // representation is (x - c)' A (x - c) <= 1, center is assumed to be origin for now MT A; @@ -66,7 +60,7 @@ class Ellipsoid{ } - // Constructor for copula ellipsoid + // Constructor for copula ellipsoid only Ellipsoid(std::vector >& Ain) { _dim = Ain.size(); A.resize(_dim, _dim); @@ -75,20 +69,6 @@ class Ellipsoid{ A(i,j) = Ain[i][j]; } } - - Eigen::SelfAdjointEigenSolver eigensolver(A); - if (eigensolver.info() != Eigen::Success) { - throw std::runtime_error("Eigen solver returned error!"); - } - - _eigen_values = eigensolver.eigenvalues(); - _Eigen_Vectors = eigensolver.eigenvectors(); - - _eigen_values_inv = _eigen_values.array().inverse().matrix(); - _eigen_values_inv_sqrt = _eigen_values_inv.array().sqrt().matrix(); - - _dim = A.rows(); - c = Point(_dim); } diff --git a/include/misc/misc.h b/include/misc/misc.h index f0a585af3..c870b871c 100644 --- a/include/misc/misc.h +++ b/include/misc/misc.h @@ -1,6 +1,7 @@ // VolEsti (volume computation and sampling library) // Copyright (c) 2012-2018 Vissarion Fisikopoulos +// Copyright (c) 2021 Vaibhav Thakkar // Licensed under GNU LGPL.3, see LICENCE file @@ -211,13 +212,10 @@ std::pair read_inner_ball(std::istream &is) { - Next `m` lines follow containing a pair 'i j' in each line to signify A_i <= A_j i.e i_th element is less than or equal to the j_th element */ -Poset read_poset_from_file(std::string filename) { +Poset read_poset_from_file(std::istream &data_file) { typedef typename Poset::RT RT; typedef typename Poset::RV RV; - std::ifstream data_file; - data_file.open(filename); - // read number of elements unsigned int n; data_file >> n; @@ -233,11 +231,9 @@ Poset read_poset_from_file(std::string filename) { // read a poset given as an adjacency matrix -std::pair read_poset_from_file_adj_matrix(std::string filename) { +std::pair read_poset_from_file_adj_matrix(std::istream &in) { typedef typename Poset::RV RV; - std::ifstream in; - in.open(filename); RV edges; unsigned int x, n = 0; From ef7ac5548542b90f30e95c69febd4bc78e28593e Mon Sep 17 00:00:00 2001 From: Vaibhav Thakkar Date: Tue, 14 Sep 2021 04:04:38 +0530 Subject: [PATCH 6/7] minor fix --- include/convex_bodies/orderpolytope.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/convex_bodies/orderpolytope.h b/include/convex_bodies/orderpolytope.h index 33f5c0d6d..b0fb8112f 100644 --- a/include/convex_bodies/orderpolytope.h +++ b/include/convex_bodies/orderpolytope.h @@ -106,7 +106,7 @@ class OrderPolytope { // print polytope in Ax <= b format - void print() + void print() const { std::cout << " " << _A.rows() << " " << _d << " double" << std::endl; for (unsigned int i = 0; i < _A.rows(); i++) { @@ -664,7 +664,7 @@ class OrderPolytope { VT inner_point() { // get topologically sorted list of indices - std::vector sorted_list = poset.topologically_sorted_list(); + std::vector sorted_list = _poset.topologically_sorted_list(); // vector to hold n linearly spaced values between 0-1 std::vector lin_space_values(_d); From e836751bb178680f9a32341e48a4b0a6337ef872 Mon Sep 17 00:00:00 2001 From: Vaibhav Thakkar Date: Tue, 14 Sep 2021 13:50:12 +0530 Subject: [PATCH 7/7] add requested changes --- include/convex_bodies/ellipsoid.h | 2 +- include/volume/{math_helpers.h => math_helpers.hpp} | 8 +++++--- include/volume/volume_cooling_balls.hpp | 2 +- include/volume/volume_cooling_gaussians.hpp | 2 +- 4 files changed, 8 insertions(+), 6 deletions(-) rename include/volume/{math_helpers.h => math_helpers.hpp} (85%) diff --git a/include/convex_bodies/ellipsoid.h b/include/convex_bodies/ellipsoid.h index 155e6db1f..a2a8a00d2 100644 --- a/include/convex_bodies/ellipsoid.h +++ b/include/convex_bodies/ellipsoid.h @@ -16,7 +16,7 @@ #include #include -#include "volume/math_helpers.h" +#include "volume/math_helpers.hpp" template diff --git a/include/volume/math_helpers.h b/include/volume/math_helpers.hpp similarity index 85% rename from include/volume/math_helpers.h rename to include/volume/math_helpers.hpp index fe1bf8155..11e2742b9 100644 --- a/include/volume/math_helpers.h +++ b/include/volume/math_helpers.hpp @@ -1,13 +1,15 @@ // VolEsti (volume computation and sampling library) +// Copyright (c) 2012-2020 Vissarion Fisikopoulos +// Copyright (c) 2018-2020 Apostolos Chalkis // Copyright (c) 2021 Vaibhav Thakkar // Contributed and/or modified by Vaibhav Thakkar, as part of Google Summer of Code 2021 program. // Licensed under GNU LGPL.3, see LICENCE file -#ifndef MATH_HELPERS_H -#define MATH_HELPERS_H +#ifndef MATH_HELPERS_HPP +#define MATH_HELPERS_HPP #include #include @@ -41,4 +43,4 @@ static NT log_gamma_function(NT x) return (std::log(x - NT(1)) + log_gamma_function(x - NT(1))); } -#endif // MATH_HELPERS_H \ No newline at end of file +#endif // MATH_HELPERS_HPP \ No newline at end of file diff --git a/include/volume/volume_cooling_balls.hpp b/include/volume/volume_cooling_balls.hpp index 5958f6f98..fda05f25b 100644 --- a/include/volume/volume_cooling_balls.hpp +++ b/include/volume/volume_cooling_balls.hpp @@ -24,7 +24,7 @@ #include "convex_bodies/ball.h" #include "convex_bodies/ballintersectconvex.h" #include "sampling/random_point_generators.hpp" -#include "math_helpers.h" +#include "volume/math_helpers.hpp" //////////////////////////////////// diff --git a/include/volume/volume_cooling_gaussians.hpp b/include/volume/volume_cooling_gaussians.hpp index aa1e2dd4c..9ec268445 100644 --- a/include/volume/volume_cooling_gaussians.hpp +++ b/include/volume/volume_cooling_gaussians.hpp @@ -23,7 +23,7 @@ #include "random_walks/gaussian_ball_walk.hpp" #include "random_walks/gaussian_cdhr_walk.hpp" #include "sampling/random_point_generators.hpp" -#include "math_helpers.h" +#include "volume/math_helpers.hpp" /////////////////// Helpers for random walks