From 7e26754bb4365420eb37650555c6d2bc3c39ba12 Mon Sep 17 00:00:00 2001 From: Panagiotis Repouskos Date: Fri, 24 May 2019 14:04:48 +0300 Subject: [PATCH 01/21] Changed class point to store coefficients to Eigen vector instead of std::vector. Made required code updates in polytopes.h, rounding.h,ballintersectconvex.h and updated CMakeLists.txt to include the directory with Eigen. --- include/cartesian_geom/point.h | 137 ++++++++++---------- include/convex_bodies/ballintersectconvex.h | 28 ++-- include/convex_bodies/polytopes.h | 21 ++- include/rounding.h | 13 +- test/CMakeLists.txt | 2 +- 5 files changed, 99 insertions(+), 102 deletions(-) diff --git a/include/cartesian_geom/point.h b/include/cartesian_geom/point.h index 0d207733c..1355f8166 100644 --- a/include/cartesian_geom/point.h +++ b/include/cartesian_geom/point.h @@ -16,85 +16,87 @@ class point { private: unsigned int d; - typedef std::vector Coeff; + + typedef Eigen::Matrix Coeff; + Coeff coeffs; typedef typename std::vector::iterator iter; public: typedef typename K::FT FT; point() {} - + point(const unsigned int dim) { d = dim; - coeffs = Coeff(d,0); + coeffs.setZero(d); } - + point(const unsigned int dim, iter begin, iter endit) { d = dim; - coeffs = Coeff(begin,endit); + coeffs.resize(d); + int i = 0; + + for (iter it=begin ; it != endit ; it++) + coeffs(i++) = *it; } - + + const Coeff& getCoefficients() const { + return coeffs; + } + int dimension() { return d; } - + void set_dimension(const unsigned int dim) { d = dim; } - + void set_coord(const unsigned int i, FT coord) { - coeffs[i] = coord; + coeffs(i) = coord; } - + FT operator[] (const unsigned int i) { - return coeffs[i]; + return coeffs(i); } - - point operator+ (point& p) { - point temp(p.dimension()); - typename Coeff::iterator tmit = temp.iter_begin(); - typename Coeff::iterator pit = p.iter_begin(); - typename Coeff::iterator mit = coeffs.begin(); - - for (; pit < p.iter_end(); ++pit, ++mit, ++tmit) { - (*tmit) = (*pit) + (*mit); - } + point operator+ (point& p) { + point temp; + temp.d = d; + temp.coeffs = coeffs + p.getCoefficients(); return temp; } - - point operator- (point& p) { - point temp(p.dimension()); - typename Coeff::iterator tmit = temp.iter_begin(); - typename Coeff::iterator pit = p.iter_begin(); - typename Coeff::iterator mit = coeffs.begin(); + void add(Coeff coeffs) { + this->coeffs = coeffs + this->coeffs; + } - for (; pit < p.iter_end(); ++pit, ++mit, ++tmit) { - (*tmit) = (*mit) - (*pit); - } + point operator- (point& p) { + point temp; + temp.d = d; + temp.coeffs = coeffs - p.getCoefficients(); return temp; } point operator* (const FT& k) { - point temp(d, iter_begin(), iter_end()); - - typename Coeff::iterator tmit = temp.iter_begin(); - - for (; tmit < temp.iter_end(); ++tmit) { - (*tmit) *= k; - } + point temp; + temp.d = d; + temp.coeffs = coeffs * k; return temp; } + point operator/ (const FT& k) { + point temp; + temp.d = d; + temp.coeffs = coeffs / k; + return temp; + } - bool operator== (point& p) { - - typename Coeff::iterator pit = p.iter_begin(); - typename Coeff::iterator mit = coeffs.begin(); + bool operator== (point& p) {//TODO check dim? + int i=0; - for ( ; pit!=p.iter_end(); ++pit, ++mit) { - if (*mit!=*pit) return false; + for (auto x : p.getCoefficients()) { + if (x != coeffs(i++)) return false; } return true; @@ -102,49 +104,40 @@ class point FT dot(point& p){ - FT res=FT(0.0); - - typename Coeff::iterator pit=p.iter_begin(); - typename Coeff::iterator mit=coeffs.begin(); - for( ; pit diff --git a/include/convex_bodies/ballintersectconvex.h b/include/convex_bodies/ballintersectconvex.h index ecdbe11a0..4337149f8 100644 --- a/include/convex_bodies/ballintersectconvex.h +++ b/include/convex_bodies/ballintersectconvex.h @@ -17,7 +17,7 @@ struct Ball{ typedef Point BallPoint; typedef typename Point::FT NT; typedef typename std::vector::iterator viterator; - + typedef Eigen::Matrix Coeff; Ball(Point cc, NT RR) : c(cc), R(RR) {} @@ -42,18 +42,20 @@ struct Ball{ std::pair line_intersect(Point r, Point v){ - viterator rit=r.iter_begin(); - viterator vit=v.iter_begin(); - viterator cit=c.iter_begin(); + Coeff r_coeffs = r.getCoefficients(); + Coeff v_coeffs = v.getCoefficients(); + Coeff c_coeffs = c.getCoefficients(); + //Point rc = r;// - _c; - viterator rcit=r.iter_begin(); +// viterator rcit=r.iter_begin(); + NT vrc(0); NT v2(0); NT rc2(0); - for( ; cit < c.iter_end() ; ++rcit, ++cit, ++rit, ++vit){ - vrc += *vit * (*rcit); - v2 += *vit * (*vit); - rc2 += *rcit * (*rcit); + for(int i=0 ; i < c.dimension() ; ++i){ + vrc += v_coeffs(i) * c_coeffs(i); + v2 += v_coeffs(i) * v_coeffs(i); + rc2 += r_coeffs(i) * r_coeffs(i); } NT disc_sqrt = std::sqrt(std::pow(vrc,2) - v2 * (rc2 - R)); @@ -66,13 +68,13 @@ struct Ball{ int rand_coord){ //Point rc = r;// - _c; - viterator rcit=r.iter_begin(); - NT vrc = *(rcit + rand_coord); + Coeff r_coeffs = r.getCoefficients(); + NT vrc = r_coeffs(0 + rand_coord) ; //NT v2 = NT(1); NT rc2(R); - for( ; rcit < r.iter_end() ; ++rcit){ - rc2 -= *rcit * (*rcit); + for(int i=0 ; i < r.dimension() ; ++i){ + rc2 -= r_coeffs(i) * r_coeffs(i); } //NT disc_sqrt = std::sqrt(std::pow(vrc,2) - v2 * (rc2 - _R)); diff --git a/include/convex_bodies/polytopes.h b/include/convex_bodies/polytopes.h index 7c048cc8d..7edcb80eb 100644 --- a/include/convex_bodies/polytopes.h +++ b/include/convex_bodies/polytopes.h @@ -297,17 +297,17 @@ class HPolytope{ //unsigned int i, j; unsigned int j; int m = num_of_hyperplanes(); - viterator rit, vit; + for (int i = 0; i < m; i++) { sum_nom = b(i); sum_denom = NT(0); j = 0; - rit = r.iter_begin(); - vit = v.iter_begin(); - for ( ; rit != r.iter_end(); rit++, vit++, j++){ - sum_nom -= A(i, j) * (*rit); - sum_denom += A(i, j) * (*vit); + VT r_coeffs = r.getCoefficients(); + VT v_coeffs = v.getCoefficients(); + for (int i=0 ; i< r.dimension(); i++, j++){ + sum_nom -= A(i, j) * r_coeffs(i); + sum_denom += A(i, j) * v_coeffs(i); } if (sum_denom == NT(0)) { //std::cout<<"div0"< rounding_min_ellipsoid(Polytope &P , std::pair Inne typedef typename Polytope::MT MT; typedef typename Polytope::VT VT; typedef typename Parameters::RNGType RNGType; + typedef Eigen::Matrix Coeff; + unsigned int n=var.n, walk_len=var.walk_steps, i, j = 0; Point c = InnerBall.first; NT radius = InnerBall.second; @@ -53,11 +55,12 @@ std::pair rounding_min_ellipsoid(Polytope &P , std::pair Inne // Store points in a matrix to call Khachiyan algorithm for the minimum volume enclosing ellipsoid boost::numeric::ublas::matrix Ap(n,randPoints.size()); typename std::list::iterator rpit=randPoints.begin(); - typename std::vector::iterator qit; + + for ( ; rpit!=randPoints.end(); rpit++, j++) { - qit = (*rpit).iter_begin(); i=0; - for ( ; qit!=(*rpit).iter_end(); qit++, i++){ - Ap(i,j)=double(*qit); + Coeff coeffs = rpit->getCoefficients(); + for (i=0 ; idimension(); i++){ + Ap(i,j)=double(coeffs(i)); } } boost::numeric::ublas::matrix Q(n,n); diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index e7a57ac1f..b491b446c 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -27,7 +27,7 @@ else () set(CMAKE_EXPORT_COMPILE_COMMANDS "ON") - #include_directories (BEFORE ../external/Eigen) + include_directories (BEFORE ../external/Eigen) include_directories (BEFORE ../external) include_directories (BEFORE ../external/minimum_ellipsoid) #include_directories (BEFORE ../include/cartesian_geom) From d945a40219af117771640b440baa4c2c4fa3b109 Mon Sep 17 00:00:00 2001 From: Panagiotis Repouskos Date: Sat, 25 May 2019 22:49:01 +0300 Subject: [PATCH 02/21] fixed bug --- include/convex_bodies/polytopes.h | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/include/convex_bodies/polytopes.h b/include/convex_bodies/polytopes.h index 7edcb80eb..6b78253c5 100644 --- a/include/convex_bodies/polytopes.h +++ b/include/convex_bodies/polytopes.h @@ -302,12 +302,11 @@ class HPolytope{ for (int i = 0; i < m; i++) { sum_nom = b(i); sum_denom = NT(0); - j = 0; VT r_coeffs = r.getCoefficients(); VT v_coeffs = v.getCoefficients(); - for (int i=0 ; i< r.dimension(); i++, j++){ - sum_nom -= A(i, j) * r_coeffs(i); - sum_denom += A(i, j) * v_coeffs(i); + for (j=0 ; j< r.dimension(); j++){ + sum_nom -= A(i, j) * r_coeffs(j); + sum_denom += A(i, j) * v_coeffs(j); } if (sum_denom == NT(0)) { //std::cout<<"div0"< Date: Wed, 7 Aug 2019 19:51:22 +0300 Subject: [PATCH 03/21] fix for quicker access to data When accessing the data of class Point don't copy the whole vector, but use the [] operator (changes requested during the pull request). --- include/convex_bodies/ballintersectconvex.h | 1 - include/convex_bodies/polytopes.h | 10 ++++------ include/rounding.h | 3 +-- 3 files changed, 5 insertions(+), 9 deletions(-) diff --git a/include/convex_bodies/ballintersectconvex.h b/include/convex_bodies/ballintersectconvex.h index 4337149f8..73269a1a7 100644 --- a/include/convex_bodies/ballintersectconvex.h +++ b/include/convex_bodies/ballintersectconvex.h @@ -47,7 +47,6 @@ struct Ball{ Coeff c_coeffs = c.getCoefficients(); //Point rc = r;// - _c; -// viterator rcit=r.iter_begin(); NT vrc(0); NT v2(0); diff --git a/include/convex_bodies/polytopes.h b/include/convex_bodies/polytopes.h index 6b78253c5..127985b1d 100644 --- a/include/convex_bodies/polytopes.h +++ b/include/convex_bodies/polytopes.h @@ -302,11 +302,9 @@ class HPolytope{ for (int i = 0; i < m; i++) { sum_nom = b(i); sum_denom = NT(0); - VT r_coeffs = r.getCoefficients(); - VT v_coeffs = v.getCoefficients(); for (j=0 ; j< r.dimension(); j++){ - sum_nom -= A(i, j) * r_coeffs(j); - sum_denom += A(i, j) * v_coeffs(j); + sum_nom -= A(i, j) * r[j]; + sum_denom += A(i, j) * v[j]; } if (sum_denom == NT(0)) { //std::cout<<"div0"< rounding_min_ellipsoid(Polytope &P , std::pair Inne for ( ; rpit!=randPoints.end(); rpit++, j++) { - Coeff coeffs = rpit->getCoefficients(); for (i=0 ; idimension(); i++){ - Ap(i,j)=double(coeffs(i)); + Ap(i,j)=double((*rpit)[i]); } } boost::numeric::ublas::matrix Q(n,n); From 7a7efdb7e6990f98a210e3b59f6980fcf5fc88dd Mon Sep 17 00:00:00 2001 From: Repouskos Panagiotis Date: Fri, 14 Feb 2020 00:32:45 +0200 Subject: [PATCH 04/21] Eigen clean existing code, make changes in leftover files with previous implementation --- R-proj/src/sample_points.cpp | 6 +-- include/cartesian_geom/point.h | 47 +++++++++++++++------ include/convex_bodies/ballintersectconvex.h | 15 ++----- include/convex_bodies/ellipsoids.h | 8 +--- include/convex_bodies/polytopes.h | 15 +++---- include/generators/polytope_generators.h | 14 +----- include/rounding.h | 2 +- include/solve_lp.h | 14 ++---- 8 files changed, 53 insertions(+), 68 deletions(-) diff --git a/R-proj/src/sample_points.cpp b/R-proj/src/sample_points.cpp index fc373e022..c9eff5428 100644 --- a/R-proj/src/sample_points.cpp +++ b/R-proj/src/sample_points.cpp @@ -327,12 +327,10 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P = R_NilValue Rcpp::NumericMatrix PointSet(dim,numpoints); typename std::list::iterator rpit=randPoints.begin(); - typename std::vector::iterator qit; unsigned int j = 0, i; for ( ; rpit!=randPoints.end(); rpit++, j++) { - qit = (*rpit).iter_begin(); i=0; - for ( ; qit!=(*rpit).iter_end(); qit++, i++){ - PointSet(i,j)=*qit; + for (i=0 ; idimension() ; i++){ + PointSet(i,j)= (*rpit)[i]; } } return PointSet; diff --git a/include/cartesian_geom/point.h b/include/cartesian_geom/point.h index 1355f8166..e5e0aa33f 100644 --- a/include/cartesian_geom/point.h +++ b/include/cartesian_geom/point.h @@ -10,6 +10,7 @@ #define POINT_H #include +#include "../external/Eigen/Eigen" template class point @@ -44,7 +45,7 @@ class point return coeffs; } - int dimension() { + int dimension() const { return d; } @@ -56,59 +57,79 @@ class point coeffs(i) = coord; } - FT operator[] (const unsigned int i) { + FT operator[] (const unsigned int i) const { return coeffs(i); } - point operator+ (point& p) { + point operator+ (const point& p) const { point temp; temp.d = d; temp.coeffs = coeffs + p.getCoefficients(); return temp; } - void add(Coeff coeffs) { + void operator+= (const point& p) { + coeffs += p.getCoefficients(); + } + + void operator+= (const Coeff& coeffs) { this->coeffs = coeffs + this->coeffs; } - point operator- (point& p) { + void operator= (const Coeff& coeffs) { + this->coeffs = coeffs; + } + + point operator- (const point& p) const { point temp; temp.d = d; temp.coeffs = coeffs - p.getCoefficients(); return temp; } - point operator* (const FT& k) { + point operator* (const FT k) const { point temp; temp.d = d; temp.coeffs = coeffs * k; return temp; } - point operator/ (const FT& k) { + void operator*= (const FT k) { + coeffs *= k; + } + + point operator/ (const FT k) const { point temp; temp.d = d; temp.coeffs = coeffs / k; return temp; } - bool operator== (point& p) {//TODO check dim? + void operator/= (const FT k) { + coeffs /= k; + } + + bool operator== (point& p) const { int i=0; - for (auto x : p.getCoefficients()) { - if (x != coeffs(i++)) return false; + for (i=0 ; icoeffs(i) != coeffs(i)) return false; } return true; } - FT dot(point& p){ + FT dot(const point& p) const { return coeffs.dot(p.getCoefficients()); } + FT dot(const Coeff& coeffs) const { + return this->coeffs.dot(coeffs); + } + - FT squared_length() { + FT squared_length() const { FT lsq = FT(0.0); @@ -118,7 +139,7 @@ class point return lsq; } - void print(){ + void print() const{ for(unsigned int i=0; i line_intersect(Point r, Point v){ - Coeff r_coeffs = r.getCoefficients(); - Coeff v_coeffs = v.getCoefficients(); - Coeff c_coeffs = c.getCoefficients(); //Point rc = r;// - _c; - NT vrc(0); - NT v2(0); - NT rc2(0); - for(int i=0 ; i < c.dimension() ; ++i){ - vrc += v_coeffs(i) * c_coeffs(i); - v2 += v_coeffs(i) * v_coeffs(i); - rc2 += r_coeffs(i) * r_coeffs(i); - } + NT vrc = v.dot(c); + NT v2 = v.dot(v); + NT rc2 = r.dot(r); + NT disc_sqrt = std::sqrt(std::pow(vrc,2) - v2 * (rc2 - R)); NT lamda1((NT(-1)*vrc + disc_sqrt)/v2); diff --git a/include/convex_bodies/ellipsoids.h b/include/convex_bodies/ellipsoids.h index 670d4ba61..e50071011 100644 --- a/include/convex_bodies/ellipsoids.h +++ b/include/convex_bodies/ellipsoids.h @@ -38,13 +38,7 @@ class copula_ellipsoid{ } NT mat_mult(Point p) { - VT q(dim); - unsigned int i = 0; - viterator pit = p.iter_begin(); - for ( ; pit!=p.iter_end(); ++pit, ++i){ - q(i)=(*pit); - } - return q.transpose()*G*q; + return p.getCoefficients().transpose()*G*p.getCoefficients(); } }; diff --git a/include/convex_bodies/polytopes.h b/include/convex_bodies/polytopes.h index 127985b1d..b3be32ba6 100644 --- a/include/convex_bodies/polytopes.h +++ b/include/convex_bodies/polytopes.h @@ -300,12 +300,9 @@ class HPolytope{ for (int i = 0; i < m; i++) { - sum_nom = b(i); - sum_denom = NT(0); - for (j=0 ; j< r.dimension(); j++){ - sum_nom -= A(i, j) * r[j]; - sum_denom += A(i, j) * v[j]; - } + sum_nom = b(i) - r.dot(A.row(i)); + sum_denom = v.dot(A.row(i)); + if (sum_denom == NT(0)) { //std::cout<<"div0"<::iterator pit; MT V(k, dim); unsigned int j; for (unsigned int i = 0; i < k; ++i) { p = get_direction(dim); - pit = p.iter_begin(); - j = 0; - for ( ; pit!=p.iter_end(); ++pit, ++j) { - V(i,j) = *pit; - } + V.row(i) = p.getCoefficients(); } Polytope VP; @@ -380,18 +375,13 @@ Polytope random_hpoly(unsigned int dim, unsigned int m) { RNGType rng(seed); boost::random::uniform_real_distribution<> urdist1(-10, 10); Point p(dim); - typename std::vector::iterator pit; MT A(m, dim); VT b(m); unsigned int j; for(unsigned int i=0; i(dim); - pit = p.iter_begin(); - j = 0; - for ( ; pit!=p.iter_end(); ++pit, ++j) { - A(i,j) = *pit; - } + A.row(i) = p.getCoefficients(); b(i) = 10.0; } diff --git a/include/rounding.h b/include/rounding.h index 9a0216656..d4259bd5b 100644 --- a/include/rounding.h +++ b/include/rounding.h @@ -12,7 +12,7 @@ #include "khach.h" -#include "Eigen" +#include "../external/Eigen/Eigen" // ----- ROUNDING ------ // // main rounding function diff --git a/include/solve_lp.h b/include/solve_lp.h index 4deec0b0a..59445f38b 100644 --- a/include/solve_lp.h +++ b/include/solve_lp.h @@ -854,10 +854,9 @@ Point PointInIntersection(MT V1, MT V2, Point direction, bool &empty) { set_add_rowmode(lp, FALSE); /* rowmode should be turned off again when done building the model */ // set the bounds - typename std::vector::iterator pit = direction.iter_begin(); - for(j=0; j Date: Fri, 14 Feb 2020 20:04:24 +0200 Subject: [PATCH 05/21] Optimizations - avoid creating copies - code cleanup - more coherent coding, using Eigen and eradicating std::vector --- include/cartesian_geom/point.h | 9 +++++++++ include/convex_bodies/ball.h | 7 +++---- include/convex_bodies/hpolytope.h | 3 +-- include/convex_bodies/vpolytope.h | 28 +++++++++------------------- 4 files changed, 22 insertions(+), 25 deletions(-) diff --git a/include/cartesian_geom/point.h b/include/cartesian_geom/point.h index bddd62cb0..3464e8b76 100644 --- a/include/cartesian_geom/point.h +++ b/include/cartesian_geom/point.h @@ -41,6 +41,11 @@ class point coeffs(i++) = *it; } + point(const Coeff& coeffs) { + d = coeffs.rows(); + this->coeffs = coeffs; + } + point(const unsigned int dim, std::vector cofs) { d = dim; iter it = cofs.begin(); @@ -51,6 +56,10 @@ class point } + void add(const Coeff& coeffs) { + this->coeffs += coeffs; + } + const Coeff& getCoefficients() const { return coeffs; } diff --git a/include/convex_bodies/ball.h b/include/convex_bodies/ball.h index 7c87ad979..9b3591c2b 100644 --- a/include/convex_bodies/ball.h +++ b/include/convex_bodies/ball.h @@ -111,10 +111,9 @@ struct Ball{ void compute_reflection (Point &v, const Point &p, const int &facet) { Point s = p; - s = s * (1.0 / std::sqrt(s.squared_length())); - s = ((-2.0 * v.dot(s)) * s); - v = s + v; - + s *= (1.0 / std::sqrt(s.squared_length())); + s *= (-2.0 * v.dot(s)); + v += s; } private: diff --git a/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h index 1fd8187b4..74e661a9d 100644 --- a/include/convex_bodies/hpolytope.h +++ b/include/convex_bodies/hpolytope.h @@ -254,8 +254,7 @@ class HPolytope{ NT sum; int m = A.rows(); for (int i = 0; i < m; i++) { - sum = b(i); - for (unsigned int j = 0; j < _d; j++) sum -= A(i, j) * p[j]; + sum = b(i) - A.row(i) * p.getCoefficients(); //Check if corresponding hyperplane is violated if (sum < NT(0)) return 0; diff --git a/include/convex_bodies/vpolytope.h b/include/convex_bodies/vpolytope.h index bc83308ac..b4cf1209c 100644 --- a/include/convex_bodies/vpolytope.h +++ b/include/convex_bodies/vpolytope.h @@ -163,15 +163,11 @@ class VPolytope{ } Point get_mean_of_vertices() { - std::vector vec(_d); - Point xc(_d), temp(_d); + Point xc(_d); for (int i = 0; i < num_of_vertices(); ++i) { - for (int j = 0; j < _d; ++j) vec[j] = V(i,j); - - temp = Point(_d, vec.begin(), vec.end()); - xc = xc + temp; + xc.add(V.row(i)); } - xc = xc * (1.0/NT(num_of_vertices())); + xc *= (1.0/NT(num_of_vertices())); return xc; } @@ -217,17 +213,15 @@ class VPolytope{ std::pair result; for (j=1; jgetCoefficients() - p0.getCoefficients(); } Bg = B; B = B.inverse(); for (i=0; i temp(_d,NT(0)); + typename std::vector::iterator pointIt; for (int i=0; i 1.0) a = -a; - a = a/a.norm(); + a /= a.norm(); Point s(_d, std::vector(&a[0], a.data()+a.cols()*a.rows())); From 5a3a9d697b6b3ceac13028e78c3cb6fada9e9754 Mon Sep 17 00:00:00 2001 From: Repouskos Panagiotis Date: Fri, 14 Feb 2020 20:44:54 +0200 Subject: [PATCH 06/21] Optimizations - Vectorize when possible the operations in line_intersect functions - Change the matrix in HPolytope to be RowMajor --- include/convex_bodies/hpolytope.h | 65 ++++++++++++++++++------------- 1 file changed, 39 insertions(+), 26 deletions(-) diff --git a/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h index 74e661a9d..16bd7d7f5 100644 --- a/include/convex_bodies/hpolytope.h +++ b/include/convex_bodies/hpolytope.h @@ -27,7 +27,7 @@ class HPolytope{ typedef typename std::vector::iterator viterator; //using RowMatrixXd = Eigen::Matrix; //typedef RowMatrixXd MT; - typedef Eigen::Matrix MT; + typedef Eigen::Matrix MT; typedef Eigen::Matrix VT; private: @@ -277,21 +277,22 @@ class HPolytope{ std::pair line_intersect(Point &r, Point &v) { NT lamda = 0, min_plus = NT(maxNT), max_minus = NT(minNT); - NT sum_nom, sum_denom; + VT sum_nom, sum_denom; //unsigned int i, j; unsigned int j; int m = num_of_hyperplanes(); - viterator rit, vit; + + + sum_nom = b - A * r.getCoefficients(); + sum_denom = A * v.getCoefficients(); for (int i = 0; i < m; i++) { - sum_nom = b(i) - A.row(i) * r.getCoefficients(); - sum_denom = A.row(i) * v.getCoefficients(); - if (sum_denom == NT(0)) { + if (sum_denom[i] == NT(0)) { //std::cout<<"div0"< 0) min_plus = lamda; if (lamda > max_minus && lamda < 0) max_minus = lamda; } @@ -306,24 +307,29 @@ class HPolytope{ std::vector &Av, bool pos = false) { NT lamda = 0, min_plus = NT(maxNT), max_minus = NT(minNT); - NT sum_nom, sum_denom, mult; + VT sum_nom, sum_denom; + NT mult; //unsigned int i, j; unsigned int j; int m = num_of_hyperplanes(), facet; viterator rit, vit, Ariter = Ar.begin(), Aviter = Av.begin(); + sum_nom.setZero(m); + sum_denom.setZero(m); + + sum_nom = b + A * r.getCoefficients(); + sum_denom = A * v.getCoefficients(); + for (int i = 0; i < m; i++, ++Ariter, ++Aviter) { - sum_nom = -A.row(i) * r.getCoefficients(); - sum_denom = A.row(i) * v.getCoefficients(); - (*Ariter) = -sum_nom; - (*Aviter) = sum_denom; - sum_nom += b(i); - if (sum_denom == NT(0)) { + (*Ariter) = -sum_nom[i]; + (*Aviter) = sum_denom[i]; + + if (sum_denom[i] == NT(0)) { //std::cout<<"div0"< 0) { min_plus = lamda; if (pos) facet = i; @@ -338,23 +344,26 @@ class HPolytope{ std::vector &Av, const NT &lambda_prev, bool pos = false) { NT lamda = 0, min_plus = NT(maxNT), max_minus = NT(minNT); - NT sum_nom, sum_denom, mult; + VT sum_denom; + NT sum_nom, mult; //unsigned int i, j; unsigned int j; int m = num_of_hyperplanes(), facet; viterator vit, Ariter = Ar.begin(), Aviter = Av.begin(); + sum_denom = A * v.getCoefficients(); + for (int i = 0; i < m; i++, ++Ariter, ++Aviter) { (*Ariter) += lambda_prev * (*Aviter); sum_nom = b(i) - (*Ariter); - sum_denom = A.row(i) * v.getCoefficients(); - (*Aviter) = sum_denom; - if (sum_denom == NT(0)) { + + (*Aviter) = sum_denom[i]; + if (sum_denom[i] == NT(0)) { //std::cout<<"div0"< 0) { min_plus = lamda; if (pos) facet = i; @@ -387,20 +396,24 @@ class HPolytope{ std::vector &lamdas) { NT lamda = 0, min_plus = NT(maxNT), max_minus = NT(minNT); - NT sum_nom, sum_denom; + VT sum_nom; + NT sum_denom; unsigned int j; int m = num_of_hyperplanes(); viterator rit; + sum_nom = b - A * r.getCoefficients(); + for (int i = 0; i < m; i++) { - sum_nom = b(i) - A.row(i)*r.getCoefficients(); sum_denom = A(i, rand_coord); - lamdas[i] = sum_nom; + + lamdas[i] = sum_nom[i]; + if (sum_denom == NT(0)) { //std::cout<<"div0"< 0) min_plus = lamda; if (lamda > max_minus && lamda < 0) max_minus = lamda; @@ -492,8 +505,8 @@ class HPolytope{ VT a = A.row(facet); Point s(_d, std::vector(&a[0], a.data()+a.cols()*a.rows())); - s = ((-2.0 * v.dot(s)) * s); - v = s + v; + s *= (-2.0 * v.dot(s)); + v += s; } From 09d33f37faaedd9f362bccbc1e7e3d08795f07a3 Mon Sep 17 00:00:00 2001 From: Repouskos Panagiotis Date: Sat, 15 Feb 2020 16:03:57 +0200 Subject: [PATCH 07/21] Major - change uses of std::vector to Eigen vector - Change the matrix in HPolytope NOT to be RowMajor (seemed slower, needs more testing) - enabled no debug macro for Eigen in /test/CmakeList - added 2 tests, to compute volume with rdhr and with BiW --- include/annealing/gaussian_annealing.h | 7 +- include/annealing/ratio_estimation.h | 13 +- include/convex_bodies/ball.h | 21 +- include/convex_bodies/ballintersectconvex.h | 15 +- include/convex_bodies/hpolytope.h | 127 +++++----- include/convex_bodies/vpolyintersectvpoly.h | 20 +- include/convex_bodies/vpolytope.h | 27 ++- include/convex_bodies/zonoIntersecthpoly.h | 20 +- include/convex_bodies/zpolytope.h | 23 +- include/lp_oracles/solve_lp.h | 7 +- include/samplers/gaussian_samplers.h | 17 +- include/samplers/samplers.h | 47 +++- include/volume/volume.h | 4 +- test/CMakeLists.txt | 8 +- test/volume_test_billiard.cpp | 252 ++++++++++++++++++++ test/volume_test_rdhr.cpp | 247 +++++++++++++++++++ 16 files changed, 710 insertions(+), 145 deletions(-) create mode 100644 test/volume_test_billiard.cpp create mode 100644 test/volume_test_rdhr.cpp diff --git a/include/annealing/gaussian_annealing.h b/include/annealing/gaussian_annealing.h index 01a1a0d33..eea190670 100644 --- a/include/annealing/gaussian_annealing.h +++ b/include/annealing/gaussian_annealing.h @@ -195,7 +195,10 @@ void get_annealing_schedule(Polytope &P, const NT &radius, const NT &ratio, cons Point p_prev=p; - std::vector lamdas(P.num_of_hyperplanes(), NT(0)); + typedef Eigen::Matrix VT; + VT lamdas; + lamdas.setZero(P.num_of_hyperplanes()); + while (true) { if (var.ball_walk) { @@ -206,7 +209,7 @@ void get_annealing_schedule(Polytope &P, const NT &radius, const NT &ratio, cons curr_fn = 0; curr_its = 0; - std::fill(lamdas.begin(), lamdas.end(), NT(0)); + lamdas.setConstant(NT(0)); steps = totalSteps; if (var.cdhr_walk){ diff --git a/include/annealing/ratio_estimation.h b/include/annealing/ratio_estimation.h index a5ac62e77..cd9e6c0a7 100644 --- a/include/annealing/ratio_estimation.h +++ b/include/annealing/ratio_estimation.h @@ -28,7 +28,12 @@ NT esti_ratio(PolyBall1 &Pb1, PolyBall2 &Pb2, const NT &ratio, const NT &error, bool print = var.verbose; NT min_val = std::numeric_limits::lowest(), max_val = std::numeric_limits::max(), val, lambda; size_t totCount = Ntot, countIn = Ntot * ratio; - std::vector last_W(W), lamdas(Pb1.num_of_hyperplanes()), Av(Pb1.num_of_hyperplanes()); + std::vector last_W(W); + typedef Eigen::Matrix VT; + VT lamdas, Av; + lamdas.setZero(Pb1.num_of_hyperplanes()); + Av.setZero(Pb1.num_of_hyperplanes()); + std::list randPoints; typename std::vector::iterator minmaxIt; typename std::list::iterator rpit; @@ -90,7 +95,11 @@ NT esti_ratio_interval(PolyBall1 &Pb1, PolyBall2 &Pb2, const NT &ratio, const NT int n = var.n, index = 0, iter = 1; bool print = var.verbose; - std::vector last_W(W), lamdas(Pb1.num_of_hyperplanes()), Av(Pb1.num_of_hyperplanes()); + std::vector last_W(W); + typedef Eigen::Matrix VT; + VT lamdas, Av; + Av.setZero(Pb1.num_of_hyperplanes()); + lamdas.setZero(Pb1.num_of_hyperplanes()); NT val, sum_sq=0.0, sum=0.0, lambda; size_t totCount = Ntot, countIn = Ntot * ratio; //std::cout<<"countIn = "< 0) min_plus = lamda; if (lamda > max_minus && lamda < 0) max_minus = lamda; } + + sum_nom_data++; + sum_denom_data++; } return std::pair(min_plus, max_minus); } @@ -303,72 +311,73 @@ class HPolytope{ // compute intersection points of a ray starting from r and pointing to v // with polytope discribed by A and b - std::pair line_intersect(Point &r, Point &v, std::vector &Ar, - std::vector &Av, bool pos = false) { + std::pair line_intersect(Point &r, Point &v, VT& Ar, + VT& Av, bool pos = false) { NT lamda = 0, min_plus = NT(maxNT), max_minus = NT(minNT); - VT sum_nom, sum_denom; + VT sum_nom; NT mult; //unsigned int i, j; unsigned int j; int m = num_of_hyperplanes(), facet; - viterator rit, vit, Ariter = Ar.begin(), Aviter = Av.begin(); - sum_nom.setZero(m); - sum_denom.setZero(m); + Ar = A * r.getCoefficients(); + sum_nom = b - Ar; + Av = A * v.getCoefficients();; - sum_nom = b + A * r.getCoefficients(); - sum_denom = A * v.getCoefficients(); - - for (int i = 0; i < m; i++, ++Ariter, ++Aviter) { - (*Ariter) = -sum_nom[i]; - (*Aviter) = sum_denom[i]; + NT* Av_data = Av.data(); + NT* sum_nom_data = sum_nom.data(); - if (sum_denom[i] == NT(0)) { + for (int i = 0; i < m; i++) { + if (*Av_data == NT(0)) { //std::cout<<"div0"< 0) { min_plus = lamda; if (pos) facet = i; }else if (lamda > max_minus && lamda < 0) max_minus = lamda; } + + Av_data++; + sum_nom_data++; } if (pos) return std::pair(min_plus, facet); return std::pair(min_plus, max_minus); } - std::pair line_intersect(Point &r, Point &v, std::vector &Ar, - std::vector &Av, const NT &lambda_prev, bool pos = false) { + std::pair line_intersect(Point &r, Point &v, VT& Ar, + VT& Av, const NT &lambda_prev, bool pos = false) { NT lamda = 0, min_plus = NT(maxNT), max_minus = NT(minNT); - VT sum_denom; - NT sum_nom, mult; + VT sum_nom; + NT mult; //unsigned int i, j; unsigned int j; int m = num_of_hyperplanes(), facet; - viterator vit, Ariter = Ar.begin(), Aviter = Av.begin(); - sum_denom = A * v.getCoefficients(); - - for (int i = 0; i < m; i++, ++Ariter, ++Aviter) { - (*Ariter) += lambda_prev * (*Aviter); - sum_nom = b(i) - (*Ariter); + Ar += lambda_prev*Av; + sum_nom = b - Ar; + Av = A * v.getCoefficients(); + NT* sum_nom_data = sum_nom.data(); + NT* Av_data = Av.data(); - (*Aviter) = sum_denom[i]; - if (sum_denom[i] == NT(0)) { + for (int i = 0; i < m; i++) { + if (*Av_data == NT(0)) { //std::cout<<"div0"< 0) { min_plus = lamda; if (pos) facet = i; }else if (lamda > max_minus && lamda < 0) max_minus = lamda; } + Av_data++; + sum_nom_data++; } if (pos) return std::pair(min_plus, facet); return std::pair(min_plus, max_minus); @@ -377,47 +386,49 @@ class HPolytope{ // compute intersection point of a ray starting from r and pointing to v // with polytope discribed by A and b - std::pair line_positive_intersect(Point &r, Point &v, std::vector &Ar, - std::vector &Av) { + std::pair line_positive_intersect(Point &r, Point &v, VT& Ar, + VT& Av) { return line_intersect(r, v, Ar, Av, true); } // compute intersection point of a ray starting from r and pointing to v // with polytope discribed by A and b - std::pair line_positive_intersect(Point &r, Point &v, std::vector &Ar, - std::vector &Av, const NT &lambda_prev) { + std::pair line_positive_intersect(Point &r, Point &v, VT& Ar, + VT& Av, const NT &lambda_prev) { return line_intersect(r, v, Ar, Av, lambda_prev, true); } //First coordinate ray intersecting convex polytope std::pair line_intersect_coord(Point &r, const unsigned int &rand_coord, - std::vector &lamdas) { + VT& lamdas) { NT lamda = 0, min_plus = NT(maxNT), max_minus = NT(minNT); - VT sum_nom; - NT sum_denom; + VT sum_denom; unsigned int j; int m = num_of_hyperplanes(); viterator rit; - sum_nom = b - A * r.getCoefficients(); + sum_denom = A.col(rand_coord); + lamdas = b - A * r.getCoefficients(); - for (int i = 0; i < m; i++) { - sum_denom = A(i, rand_coord); + NT* lamda_data = lamdas.data(); + NT* sum_denom_data = sum_denom.data(); - lamdas[i] = sum_nom[i]; + for (int i = 0; i < m; i++) { - if (sum_denom == NT(0)) { + if (*sum_denom_data == NT(0)) { //std::cout<<"div0"< 0) min_plus = lamda; if (lamda > max_minus && lamda < 0) max_minus = lamda; } + lamda_data++; + sum_denom_data++; } return std::pair(min_plus, max_minus); } @@ -428,29 +439,28 @@ class HPolytope{ const Point &r_prev, const unsigned int rand_coord, const unsigned int rand_coord_prev, - std::vector &lamdas) { - - viterator lamdait = lamdas.begin(); + VT& lamdas) { + ; NT lamda = 0, min_plus = NT(maxNT), max_minus = NT(minNT); - NT sum_nom, sum_denom, c_rand_coord, c_rand_coord_prev; + int m = num_of_hyperplanes(); + lamdas += A.col(rand_coord_prev)* (r_prev[rand_coord_prev] - r[rand_coord_prev]); + NT* data = lamdas.data(); + for (int i = 0; i < m; i++) { - sum_denom = b(i); - c_rand_coord = A(i, rand_coord); - c_rand_coord_prev = A(i, rand_coord_prev); + NT a = A(i, rand_coord); - *lamdait = *lamdait + c_rand_coord_prev * (r_prev[rand_coord_prev] - r[rand_coord_prev]); - if (c_rand_coord == NT(0)) { + if (a == NT(0)) { //std::cout<<"div0"< 0) min_plus = lamda; if (lamda > max_minus && lamda < 0) max_minus = lamda; } - ++lamdait; + data++; } return std::pair(min_plus, max_minus); } @@ -503,8 +513,9 @@ class HPolytope{ void compute_reflection(Point &v, const Point &p, const int facet) { - VT a = A.row(facet); - Point s(_d, std::vector(&a[0], a.data()+a.cols()*a.rows())); +// VT a = A.row(facet); +// Point s(_d, std::vector(&a[0], a.data()+a.cols()*a.rows())); + Point s(A.row(facet)); s *= (-2.0 * v.dot(s)); v += s; diff --git a/include/convex_bodies/vpolyintersectvpoly.h b/include/convex_bodies/vpolyintersectvpoly.h index cc73ac1b9..4b5b09a97 100644 --- a/include/convex_bodies/vpolyintersectvpoly.h +++ b/include/convex_bodies/vpolyintersectvpoly.h @@ -210,16 +210,16 @@ class IntersectionOfVpoly { // compute intersection point of ray starting from r and pointing to v // with the V-polytope - std::pair line_intersect(const Point &r, const Point &v, const std::vector &Ar, - const std::vector &Av) { + std::pair line_intersect(const Point &r, const Point &v, const VT &Ar, + const VT &Av) { return line_intersect(r, v); } // compute intersection point of ray starting from r and pointing to v // with the V-polytope - std::pair line_intersect(const Point &r, const Point &v, const std::vector &Ar, - const std::vector &Av, const NT &lambda) { + std::pair line_intersect(const Point &r, const Point &v, const VT &Ar, + const VT &Av, const NT &lambda) { return line_intersect(r, v); } @@ -234,14 +234,14 @@ class IntersectionOfVpoly { return std::pair(P2pair.first, 2); } - std::pair line_positive_intersect(const Point &r, const Point &v, const std::vector &Ar, - const std::vector &Av) { + std::pair line_positive_intersect(const Point &r, const Point &v, const VT &Ar, + const VT &Av) { return line_positive_intersect(r, v); } - std::pair line_positive_intersect(const Point &r, const Point &v, const std::vector &Ar, - const std::vector &Av, const NT &lambda_prev) { + std::pair line_positive_intersect(const Point &r, const Point &v, const VT &Ar, + const VT &Av, const NT &lambda_prev) { return line_positive_intersect(r, v);//, Ar, Av); } @@ -250,7 +250,7 @@ class IntersectionOfVpoly { // with the V-polytope std::pair line_intersect_coord(const Point &r, const unsigned int &rand_coord, - const std::vector &lamdas) { + const VT &lamdas) { std::pair P1pair = P1.line_intersect_coord(r, rand_coord, lamdas); std::pair P2pair = P2.line_intersect_coord(r, rand_coord, lamdas); return std::pair(std::min(P1pair.first, P2pair.first), @@ -264,7 +264,7 @@ class IntersectionOfVpoly { const Point &r_prev, const unsigned int &rand_coord, const unsigned int &rand_coord_prev, - const std::vector &lamdas) { + const VT &lamdas) { return line_intersect_coord(r, rand_coord, lamdas); } diff --git a/include/convex_bodies/vpolytope.h b/include/convex_bodies/vpolytope.h index b4cf1209c..0d2cb5393 100644 --- a/include/convex_bodies/vpolytope.h +++ b/include/convex_bodies/vpolytope.h @@ -295,10 +295,15 @@ class VPolytope{ boost::numeric::ublas::matrix Ap(_d,randPoints.size()); typename std::list::iterator rpit=randPoints.begin(); + + unsigned int i, j = 0; for ( ; rpit!=randPoints.end(); rpit++, j++) { + const NT* point_data = rpit->getCoefficients().data(); + for ( i=0; i < rpit->dimension(); i++){ - Ap(i,j)=double((*rpit)[i]); + Ap(i,j)=double(*point_data); + point_data++; } } boost::numeric::ublas::matrix Q(_d, _d); @@ -343,16 +348,16 @@ class VPolytope{ // compute intersection point of ray starting from r and pointing to v // with the V-polytope - std::pair line_intersect(const Point &r, const Point &v, const std::vector &Ar, - const std::vector &Av) { + std::pair line_intersect(const Point &r, const Point &v, const VT &Ar, + const VT &Av) { return intersect_double_line_Vpoly(V, r, v, row, colno); } // compute intersection point of ray starting from r and pointing to v // with the V-polytope - std::pair line_intersect(const Point &r, const Point &v, const std::vector &Ar, - const std::vector &Av, const NT &lambda_prev) { + std::pair line_intersect(const Point &r, const Point &v, const VT &Ar, + const VT &Av, const NT &lambda_prev) { return intersect_double_line_Vpoly(V, r, v, row, colno); } @@ -362,14 +367,14 @@ class VPolytope{ return std::pair (intersect_line_Vpoly(V, r, v, conv_comb, row, colno, false, false), 1); } - std::pair line_positive_intersect(const Point &r, const Point &v, const std::vector &Ar, - const std::vector &Av) { + std::pair line_positive_intersect(const Point &r, const Point &v, const VT &Ar, + const VT &Av) { return line_positive_intersect(r, v); } - std::pair line_positive_intersect(const Point &r, const Point &v, const std::vector &Ar, - const std::vector &Av, const NT &lambda_prev) { + std::pair line_positive_intersect(const Point &r, const Point &v, const VT &Ar, + const VT &Av, const NT &lambda_prev) { return line_positive_intersect(r, v);//, Ar, Av); } @@ -378,7 +383,7 @@ class VPolytope{ // with the V-polytope std::pair line_intersect_coord(const Point &r, const unsigned int rand_coord, - const std::vector &lamdas) { + const VT &lamdas) { std::vector temp(_d); temp[rand_coord]=1.0; @@ -393,7 +398,7 @@ class VPolytope{ const Point &r_prev, const unsigned int rand_coord, const unsigned int rand_coord_prev, - const std::vector &lamdas) { + const VT &lamdas) { return line_intersect_coord(r, rand_coord, lamdas); } diff --git a/include/convex_bodies/zonoIntersecthpoly.h b/include/convex_bodies/zonoIntersecthpoly.h index c39182efb..08ead19d1 100644 --- a/include/convex_bodies/zonoIntersecthpoly.h +++ b/include/convex_bodies/zonoIntersecthpoly.h @@ -63,16 +63,16 @@ class ZonoIntersectHPoly { std::max(polypair.second, zonopair.second)); } - std::pair line_intersect(Point &r, Point &v, std::vector &Ar, - std::vector &Av) { + std::pair line_intersect(Point &r, Point &v, VT &Ar, + VT &Av) { std::pair polypair = HP.line_intersect(r, v, Ar, Av); std::pair zonopair = Z.line_intersect(r, v); return std::pair(std::min(polypair.first, zonopair.first), std::max(polypair.second, zonopair.second)); } - std::pair line_intersect(Point &r, Point &v, std::vector &Ar, - std::vector &Av, NT &lambda_prev) { + std::pair line_intersect(Point &r, Point &v, VT &Ar, + VT &Av, NT &lambda_prev) { std::pair polypair = HP.line_intersect(r, v, Ar, Av, lambda_prev); std::pair zonopair = Z.line_intersect(r, v); return std::pair(std::min(polypair.first, zonopair.first), @@ -81,7 +81,7 @@ class ZonoIntersectHPoly { //First coordinate ray shooting intersecting convex body std::pair line_intersect_coord(Point &r,const unsigned int &rand_coord, - std::vector &lamdas) { + VT &lamdas) { std::pair polypair = HP.line_intersect_coord(r, rand_coord, lamdas); std::pair zonopair = Z.line_intersect_coord(r, rand_coord, lamdas); @@ -90,8 +90,8 @@ class ZonoIntersectHPoly { } - std::pair line_positive_intersect(Point &r, Point &v, std::vector &Ar, - std::vector &Av) { + std::pair line_positive_intersect(Point &r, Point &v, VT &Ar, + VT &Av) { std::pair polypair = HP.line_positive_intersect(r, v, Ar, Av); std::pair zonopair = Z.line_positive_intersect(r, v, Ar, Av); @@ -102,8 +102,8 @@ class ZonoIntersectHPoly { return std::pair(std::min(polypair.first, zonopair.first), facet); } - std::pair line_positive_intersect(Point &r, Point &v, std::vector &Ar, - std::vector &Av, NT &lambda_prev) { + std::pair line_positive_intersect(Point &r, Point &v, VT &Ar, + VT &Av, NT &lambda_prev) { std::pair polypair = HP.line_positive_intersect(r, v, Ar, Av, lambda_prev); std::pair zonopair = Z.line_positive_intersect(r, v, Ar, Av); int facet = HP.num_of_hyperplanes()+1; @@ -118,7 +118,7 @@ class ZonoIntersectHPoly { const Point &r_prev, const unsigned int &rand_coord, const unsigned int &rand_coord_prev, - std::vector &lamdas) { + VT &lamdas) { std::pair polypair = HP.line_intersect_coord(r, r_prev, rand_coord, rand_coord_prev, lamdas); std::pair zonopair = Z.line_intersect_coord(r, rand_coord, lamdas); diff --git a/include/convex_bodies/zpolytope.h b/include/convex_bodies/zpolytope.h index 98fef55d7..f074d58c3 100644 --- a/include/convex_bodies/zpolytope.h +++ b/include/convex_bodies/zpolytope.h @@ -288,27 +288,27 @@ class Zonotope { // compute intersection point of ray starting from r and pointing to v // with the Zonotope - std::pair line_intersect(const Point &r, const Point &v, const std::vector &Ar, - const std::vector &Av) { + std::pair line_intersect(const Point &r, const Point &v, const VT& Ar, + const VT& Av) { return intersect_line_zono(V, r, v, conv_comb, colno); } // compute intersection point of ray starting from r and pointing to v // with the Zonotope - std::pair line_intersect(const Point &r, const Point &v, const std::vector &Ar, - const std::vector &Av, const NT &lambda_prev) { + std::pair line_intersect(const Point &r, const Point &v, const VT &Ar, + const VT &Av, const NT &lambda_prev) { return intersect_line_zono(V, r, v, conv_comb, colno); } - std::pair line_positive_intersect(const Point &r, const Point &v, const std::vector &Ar, - const std::vector &Av) { + std::pair line_positive_intersect(const Point &r, const Point &v, const VT &Ar, + const VT &Av) { return std::pair (intersect_line_Vpoly(V, r, v, conv_comb, row, colno, false, true), 1); } - std::pair line_positive_intersect(const Point &r, const Point &v, const std::vector &Ar, - const std::vector &Av, const NT &lambda_prev) { + std::pair line_positive_intersect(const Point &r, const Point &v, const VT &Ar, + const VT &Av, const NT &lambda_prev) { return line_positive_intersect(r, v, Ar, Av); } @@ -317,7 +317,7 @@ class Zonotope { // with the Zonotope std::pair line_intersect_coord(const Point &r, const unsigned int rand_coord, - const std::vector &lamdas) { + const VT &lamdas) { std::vector temp(_d,0); temp[rand_coord]=1.0; @@ -334,7 +334,7 @@ class Zonotope { const Point &r_prev, const unsigned int rand_coord, const unsigned int rand_coord_prev, - const std::vector &lamdas) { + const VT &lamdas) { return line_intersect_coord(r, rand_coord, lamdas); } @@ -383,8 +383,7 @@ class Zonotope { } VT a = Fmat.fullPivLu().kernel(); - NT sum = 0.0; - for (int k = 0; k < _d; ++k) sum += a(k)*p[k]; + NT sum = p.getCoefficients().dot(a); if(sum<0.0) a = -1.0*a; diff --git a/include/lp_oracles/solve_lp.h b/include/lp_oracles/solve_lp.h index 310f578ca..a3224e67f 100644 --- a/include/lp_oracles/solve_lp.h +++ b/include/lp_oracles/solve_lp.h @@ -242,10 +242,15 @@ Point PointInIntersection(MT V1, MT V2, Point direction, bool &empty) { set_add_rowmode(lp, FALSE); /* rowmode should be turned off again when done building the model */ + const NT* direction_data = direction.getCoefficients().data(); + REAL* row_temp = row; + // set the bounds for(j=0; j +template void gaussian_first_coord_point(Polytope &P, Point &p, // a point to start Point &p_prev, // previous point unsigned int &coord_prev, // previous coordinate ray unsigned int walk_len, // number of steps for the random walk const NT &a_i, - std::vector &lamdas, + VT &lamdas, Parameters const& var) { typedef typename Parameters::RNGType RNGType; unsigned int n = var.n, rand_coord; @@ -175,14 +175,14 @@ void gaussian_first_coord_point(Polytope &P, // Compute the next point with target distribution the gaussian -template +template void gaussian_next_point(Polytope &P, Point &p, // a point to start Point &p_prev, // previous point unsigned int &coord_prev, // previous coordinate ray const unsigned int walk_len, // number of steps for the random walk const NT &a_i, - std::vector &lamdas, + VT &lamdas, Parameters const& var) { typedef typename Parameters::RNGType RNGType; unsigned int n = var.n, rand_coord; @@ -223,7 +223,10 @@ void rand_gaussian_point_generator(Polytope &P, RNGType &rng2 = var.rng; boost::random::uniform_int_distribution<> uidist(0, n - 1); - std::vector lamdas(P.num_of_hyperplanes(), NT(0)); + typedef Eigen::Matrix VT; + VT lamdas; + lamdas.setZero(P.num_of_hyperplanes()); + unsigned int rand_coord = uidist(rng2), coord_prev, rand_coord_prev; NT ball_rad = var.delta; Point p_prev = p; @@ -282,14 +285,14 @@ void gaussian_hit_and_run(Point &p, // hit-and-run with orthogonal directions and update -template +template void gaussian_hit_and_run_coord_update(Point &p, Point &p_prev, Polytope &P, unsigned int rand_coord, unsigned int rand_coord_prev, const NT &a_i, - std::vector &lamdas, + VT &lamdas, Parameters const& var) { std::pair bpair = P.line_intersect_coord(p, p_prev, rand_coord, rand_coord_prev, lamdas); NT dis = rand_exp_range_coord(p[rand_coord] + bpair.second, p[rand_coord] + bpair.first, a_i, var); diff --git a/include/samplers/samplers.h b/include/samplers/samplers.h index 88462c3c0..c479eb9ea 100644 --- a/include/samplers/samplers.h +++ b/include/samplers/samplers.h @@ -11,6 +11,8 @@ #define RANDOM_SAMPLERS_H + + // Pick a random direction as a normilized vector template Point get_direction(const unsigned int dim) { @@ -111,12 +113,20 @@ void boundary_rand_point_generator(Polytope &P, { typedef typename Parameters::RNGType RNGType; typedef typename Point::FT NT; + typedef Eigen::Matrix VT; unsigned int n = var.n; RNGType &rng = var.rng; boost::random::uniform_real_distribution<> urdist(0, 1); boost::random::uniform_int_distribution<> uidist(0, n - 1); - std::vector lamdas(P.num_of_hyperplanes(), NT(0)), Av(P.num_of_hyperplanes(), NT(0)); + typedef Eigen::Matrix VT; + VT lamdas, Av; + + lamdas.setZero(P.num_of_hyperplanes()); + Av.setZero(P.num_of_hyperplanes()); + + + unsigned int rand_coord, rand_coord_prev; NT kapa, lambda; Point p_prev = p, p1(n), p2(n), v(n); @@ -184,7 +194,12 @@ void rand_point_generator(Polytope &P, boost::random::uniform_real_distribution<> urdist(0, 1); boost::random::uniform_int_distribution<> uidist(0, n - 1); - std::vector lamdas(P.num_of_hyperplanes(), NT(0)), Av(P.num_of_hyperplanes(), NT(0)); + typedef Eigen::Matrix VT; + VT lamdas, Av; + + lamdas.setZero(P.num_of_hyperplanes()); + Av.setZero(P.num_of_hyperplanes()); + unsigned int rand_coord, rand_coord_prev; NT kapa, ball_rad = var.delta, lambda; Point p_prev = p, v(n); @@ -249,7 +264,13 @@ void rand_point_generator(BallPoly &PBLarge, RNGType &rng = var.rng; boost::random::uniform_real_distribution<> urdist(0, 1); boost::random::uniform_int_distribution<> uidist(0, n - 1); - std::vector lamdas(PBLarge.num_of_hyperplanes(), NT(0)), Av(PBLarge.num_of_hyperplanes(), NT(0)); + + typedef Eigen::Matrix VT; + VT lamdas, Av; + + lamdas.setZero(PBLarge.num_of_hyperplanes()); + Av.setZero(PBLarge.num_of_hyperplanes()); + unsigned int rand_coord, rand_coord_prev; NT kapa, ball_rad = var.delta, lambda; Point p_prev = p, v(n); @@ -299,14 +320,14 @@ void rand_point_generator(BallPoly &PBLarge, } } -template +template void uniform_first_point(Polytope &P, Point &p, // a point to start Point &p_prev, // previous point unsigned int &coord_prev, // previous coordinate ray unsigned int walk_len, // number of steps for the random walk - std::vector &lamdas, - std::vector &Av, + VT &lamdas, + VT &Av, NT &lambda, const Parameters &var) { typedef typename Parameters::RNGType RNGType; @@ -363,14 +384,14 @@ void uniform_first_point(Polytope &P, -template +template void uniform_next_point(Polytope &P, Point &p, // a point to start Point &p_prev, // previous point unsigned int &coord_prev, // previous coordinate ray const unsigned int walk_len, // number of steps for the random walk - std::vector &lamdas, - std::vector &Av, + VT &lamdas, + VT &Av, NT &lambda, const Parameters &var) { typedef typename Parameters::RNGType RNGType; @@ -426,22 +447,22 @@ void hit_and_run(Point &p, //hit-and-run with orthogonal directions and update -template +template void hit_and_run_coord_update(Point &p, Point &p_prev, Polytope &P, unsigned int rand_coord, unsigned int rand_coord_prev, const NT &kapa, - std::vector &lamdas) { + VT &lamdas) { std::pair bpair = P.line_intersect_coord(p, p_prev, rand_coord, rand_coord_prev, lamdas); p_prev = p; p.set_coord(rand_coord, p[rand_coord] + bpair.first + kapa * (bpair.second - bpair.first)); } -template -void billiard_walk(ConvexBody &P, Point &p, NT diameter, std::vector &Ar, std::vector &Av, NT &lambda_prev, +template +void billiard_walk(ConvexBody &P, Point &p, NT diameter, VT &Ar, VT &Av, NT &lambda_prev, Parameters &var, bool first = false) { typedef typename Parameters::RNGType RNGType; diff --git a/include/volume/volume.h b/include/volume/volume.h index c7dac676e..58923fb6e 100644 --- a/include/volume/volume.h +++ b/include/volume/volume.h @@ -319,7 +319,9 @@ NT volume_gaussian_annealing(Polytope &P, // Initialization for the approximation of the ratios unsigned int W = var.W, coord_prev, i=0; - std::vector last_W2(W,0), fn(mm,0), its(mm,0), lamdas(m,0); + std::vector last_W2(W,0), fn(mm,0), its(mm,0); + VT lamdas; + lamdas.setZero(m); vol=std::pow(M_PI/a_vals[0], (NT(n))/2.0)*std::abs(round_value); Point p(n), p_prev(n); // The origin is the Chebychev center of the Polytope viterator fnIt = fn.begin(), itsIt = its.begin(), avalsIt = a_vals.begin(), minmaxIt; diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 05e7d1c29..0e115eace 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -50,7 +50,8 @@ else () include_directories (BEFORE ../include/lp_oracles) include_directories (BEFORE ../include/misc) - + #for Eigen + add_compile_definitions("EIGEN_NO_DEBUG") add_definitions(${CMAKE_CXX_FLAGS} "-std=c++11") # enable C++11 standard add_definitions(${CMAKE_CXX_FLAGS} "-O3") # optimization of the compiler #add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lgsl") @@ -67,6 +68,9 @@ else () add_library(test_main OBJECT test_main.cpp) add_executable (volume_test volume_test.cpp $) + add_executable (volume_test_billiard volume_test_billiard.cpp $) + add_executable (volume_test_rdhr volume_test_rdhr.cpp $) + add_executable (cheb_test chebychev_test.cpp $) #add_executable (rounding_test rounding_test.cpp $) add_executable (volumeCV_test volumeCV_test.cpp $) @@ -115,6 +119,8 @@ else () #TARGET_LINK_LIBRARIES(volume ${LP_SOLVE}) TARGET_LINK_LIBRARIES(generate ${LP_SOLVE}) TARGET_LINK_LIBRARIES(volume_test ${LP_SOLVE}) + TARGET_LINK_LIBRARIES(volume_test_rdhr ${LP_SOLVE}) + TARGET_LINK_LIBRARIES(volume_test_billiard ${LP_SOLVE}) TARGET_LINK_LIBRARIES(cheb_test ${LP_SOLVE}) #TARGET_LINK_LIBRARIES(rounding_test ${LP_SOLVE}) TARGET_LINK_LIBRARIES(volumeCV_test ${LP_SOLVE}) diff --git a/test/volume_test_billiard.cpp b/test/volume_test_billiard.cpp new file mode 100644 index 000000000..1a4f7570c --- /dev/null +++ b/test/volume_test_billiard.cpp @@ -0,0 +1,252 @@ +// VolEsti (volume computation and sampling library) + +// Copyright (c) 20012-2018 Vissarion Fisikopoulos +// Copyright (c) 2018 Apostolos Chalkis + +// Licensed under GNU LGPL.3, see LICENCE file + +#include "doctest.h" +#include +#include "Eigen/Eigen" +#include +#include "random.hpp" +#include "random/uniform_int.hpp" +#include "random/normal_distribution.hpp" +#include "random/uniform_real_distribution.hpp" +#include "volume.h" +#include "misc.h" +#include "polytope_generators.h" +#include + +template +NT factorial(NT n) +{ + return (n == 1 || n == 0) ? 1 : factorial(n - 1) * n; +} + +template +void test_volume(Polytope &HP, NT expected, NT tolerance=0.1) +{ + + typedef typename Polytope::PolytopePoint Point; + + // Setup the parameters + int n = HP.dimension(); + int walk_len=10 + n/10; + int nexp=1, n_threads=1; + NT e=1, err=0.0000000001; + int rnum = std::pow(e,-2) * 400 * n * std::log(n); + unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); + RNGType rng(seed); + boost::normal_distribution<> rdist(0,1); + boost::random::uniform_real_distribution<>(urdist); + boost::random::uniform_real_distribution<> urdist1(-1,1); + + std::pair InnerBall; + InnerBall = HP.ComputeInnerBall(); + NT diameter; + HP.comp_diam(diameter, InnerBall.second); + + vars var(rnum,n,walk_len,n_threads,err,e,0,0,0,InnerBall.second,diameter,rng, + urdist,urdist1,-1.0,false,false,false,false,false,false,false,false,true); + + //Compute chebychev ball// + std::pair CheBall; + + // Estimate the volume + std::cout << "Number type: " << typeid(NT).name() << std::endl; + NT vol = 0; + unsigned int const num_of_exp = 10; + for (unsigned int i=0; i(10, false); + test_volume(P, 1024.0); + + std::cout << "--- Testing volume of H-cube20" << std::endl; + P = gen_cube(20, false); + test_volume(P, 1048576.0); + + std::cout << "--- Testing volume of H-cube30" << std::endl; + P = gen_cube(30, false); + test_volume(P, 1073742000.0, 0.2); +} + +template +void call_test_cross(){ + typedef Cartesian Kernel; + typedef typename Kernel::Point Point; + typedef boost::mt19937 RNGType; + typedef HPolytope Hpolytope; + + std::cout << "--- Testing volume of H-cross10" << std::endl; + Hpolytope P = gen_cross(10, false); + test_volume(P, 0.0002821869); +} + +template +void call_test_birk() { + typedef Cartesian Kernel; + typedef typename Kernel::Point Point; + typedef boost::mt19937 RNGType; + typedef HPolytope Hpolytope; + Hpolytope P; + + std::cout << "--- Testing volume of H-birk3" << std::endl; + std::ifstream inp; + std::vector > Pin; + inp.open("../R-proj/inst/extdata/birk3.ine",std::ifstream::in); + read_pointset(inp,Pin); + P.init(Pin); + test_volume(P, 0.125); + + std::cout << "--- Testing volume of H-birk4" << std::endl; + std::ifstream inp2; + std::vector > Pin2; + inp2.open("../R-proj/inst/extdata/birk4.ine",std::ifstream::in); + read_pointset(inp2,Pin2); + P.init(Pin2); + test_volume(P, 0.000970018); + + std::cout << "--- Testing volume of H-birk5" << std::endl; + std::ifstream inp3; + std::vector > Pin3; + inp3.open("../R-proj/inst/extdata/birk5.ine",std::ifstream::in); + read_pointset(inp3,Pin3); + P.init(Pin3); + test_volume(P, 0.000000225); + + std::cout << "--- Testing volume of H-birk6" << std::endl; + std::ifstream inp4; + std::vector > Pin4; + inp4.open("../R-proj/inst/extdata/birk6.ine",std::ifstream::in); + read_pointset(inp4,Pin4); + P.init(Pin4); + test_volume(P, 0.0000000000009455459196, 0.5); +} + +template +void call_test_prod_simplex() { + typedef Cartesian Kernel; + typedef typename Kernel::Point Point; + typedef boost::mt19937 RNGType; + typedef HPolytope Hpolytope; + Hpolytope P; + + std::cout << "--- Testing volume of H-prod_simplex5" << std::endl; + P = gen_prod_simplex(5); + test_volume(P, std::pow(1.0 / factorial(5.0), 2)); + + std::cout << "--- Testing volume of H-prod_simplex10" << std::endl; + P = gen_prod_simplex(10); + test_volume(P, std::pow(1.0 / factorial(10.0), 2)); + + std::cout << "--- Testing volume of H-prod_simplex15" << std::endl; + P = gen_prod_simplex(15); + test_volume(P, std::pow(1.0 / factorial(15.0), 2)); + + std::cout << "--- Testing volume of H-prod_simplex20" << std::endl; + P = gen_prod_simplex(20); + test_volume(P, std::pow(1.0 / factorial(20.0), 2)); +} + +template +void call_test_simplex() { + typedef Cartesian Kernel; + typedef typename Kernel::Point Point; + typedef boost::mt19937 RNGType; + typedef HPolytope Hpolytope; + Hpolytope P; + + std::cout << "--- Testing volume of H-simplex10" << std::endl; + P = gen_simplex(10, false); + test_volume(P, 1.0 / factorial(10.0)); + + std::cout << "--- Testing volume of H-simplex20" << std::endl; + P = gen_simplex(20, false); + test_volume(P, 1.0 / factorial(20.0)); + + std::cout << "--- Testing volume of H-simplex30" << std::endl; + P = gen_simplex(30, false); + test_volume(P, 1.0 / factorial(30.0)); + + std::cout << "--- Testing volume of H-simplex40" << std::endl; + P = gen_simplex(40, false); + test_volume(P, 1.0 / factorial(40.0)); + + //std::cout << "--- Testing volume of H-simplex50" << std::endl; + //P = gen_simplex(50, false); + //test_volume(P, 1.0 / factorial(50.0)); +} + +template +void call_test_skinny_cube() { + typedef Cartesian Kernel; + typedef typename Kernel::Point Point; + typedef boost::mt19937 RNGType; + typedef HPolytope Hpolytope; + Hpolytope P; + + std::cout << "--- Testing volume of H-skinny_cube10" << std::endl; + P = gen_skinny_cube(10); + test_volume(P, 102400.0); + + //std::cout << "--- Testing volume of H-skinny_cube20" << std::endl; + //P = gen_skinny_cube(20); + //test_volume(P, 104857600.0); +} + + +TEST_CASE("cube") { + call_test_cube(); + //call_test_cube(); + //call_test_cube(); +} + +TEST_CASE("cross") { + call_test_cross(); + //call_test_cross(); + //call_test_cross(); +} + +TEST_CASE("birk") { + call_test_birk(); + //call_test_birk(); + //call_test_birk(); +} + +TEST_CASE("prod_simplex") { + call_test_prod_simplex(); + //call_test_prod_simplex(); + //call_test_prod_simplex(); +} + +TEST_CASE("simplex") { + call_test_simplex(); + //call_test_simplex(); + //call_test_simplex(); +} + +TEST_CASE("skinny_cube") { + call_test_skinny_cube(); + //call_test_skinny_cube(); + //call_test_skinny_cube(); +} diff --git a/test/volume_test_rdhr.cpp b/test/volume_test_rdhr.cpp new file mode 100644 index 000000000..15ef868e3 --- /dev/null +++ b/test/volume_test_rdhr.cpp @@ -0,0 +1,247 @@ +// VolEsti (volume computation and sampling library) + +// Copyright (c) 20012-2018 Vissarion Fisikopoulos +// Copyright (c) 2018 Apostolos Chalkis + +// Licensed under GNU LGPL.3, see LICENCE file + +#include "doctest.h" +#include +#include "Eigen/Eigen" +#include +#include "random.hpp" +#include "random/uniform_int.hpp" +#include "random/normal_distribution.hpp" +#include "random/uniform_real_distribution.hpp" +#include "volume.h" +#include "misc.h" +#include "polytope_generators.h" +#include + +template +NT factorial(NT n) +{ + return (n == 1 || n == 0) ? 1 : factorial(n - 1) * n; +} + +template +void test_volume(Polytope &HP, NT expected, NT tolerance=0.1) +{ + + typedef typename Polytope::PolytopePoint Point; + + // Setup the parameters + int n = HP.dimension(); + int walk_len=10 + n/10; + int nexp=1, n_threads=1; + NT e=1, err=0.0000000001; + int rnum = std::pow(e,-2) * 400 * n * std::log(n); + unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); + RNGType rng(seed); + boost::normal_distribution<> rdist(0,1); + boost::random::uniform_real_distribution<>(urdist); + boost::random::uniform_real_distribution<> urdist1(-1,1); + + vars var(rnum,n,walk_len,n_threads,err,e,0,0,0,0,0.0,rng, + urdist,urdist1,-1.0,false,false,false,false,false,false,true, false,false); + + //Compute chebychev ball// + std::pair CheBall; + + // Estimate the volume + std::cout << "Number type: " << typeid(NT).name() << std::endl; + NT vol = 0; + unsigned int const num_of_exp = 10; + for (unsigned int i=0; i(10, false); + test_volume(P, 1024.0); + + std::cout << "--- Testing volume of H-cube20" << std::endl; + P = gen_cube(20, false); + test_volume(P, 1048576.0); + + std::cout << "--- Testing volume of H-cube30" << std::endl; + P = gen_cube(30, false); + test_volume(P, 1073742000.0, 0.2); +} + +template +void call_test_cross(){ + typedef Cartesian Kernel; + typedef typename Kernel::Point Point; + typedef boost::mt19937 RNGType; + typedef HPolytope Hpolytope; + + std::cout << "--- Testing volume of H-cross10" << std::endl; + Hpolytope P = gen_cross(10, false); + test_volume(P, 0.0002821869); +} + +template +void call_test_birk() { + typedef Cartesian Kernel; + typedef typename Kernel::Point Point; + typedef boost::mt19937 RNGType; + typedef HPolytope Hpolytope; + Hpolytope P; + + std::cout << "--- Testing volume of H-birk3" << std::endl; + std::ifstream inp; + std::vector > Pin; + inp.open("../R-proj/inst/extdata/birk3.ine",std::ifstream::in); + read_pointset(inp,Pin); + P.init(Pin); + test_volume(P, 0.125); + + std::cout << "--- Testing volume of H-birk4" << std::endl; + std::ifstream inp2; + std::vector > Pin2; + inp2.open("../R-proj/inst/extdata/birk4.ine",std::ifstream::in); + read_pointset(inp2,Pin2); + P.init(Pin2); + test_volume(P, 0.000970018); + + std::cout << "--- Testing volume of H-birk5" << std::endl; + std::ifstream inp3; + std::vector > Pin3; + inp3.open("../R-proj/inst/extdata/birk5.ine",std::ifstream::in); + read_pointset(inp3,Pin3); + P.init(Pin3); + test_volume(P, 0.000000225); + + std::cout << "--- Testing volume of H-birk6" << std::endl; + std::ifstream inp4; + std::vector > Pin4; + inp4.open("../R-proj/inst/extdata/birk6.ine",std::ifstream::in); + read_pointset(inp4,Pin4); + P.init(Pin4); + test_volume(P, 0.0000000000009455459196, 0.5); +} + +template +void call_test_prod_simplex() { + typedef Cartesian Kernel; + typedef typename Kernel::Point Point; + typedef boost::mt19937 RNGType; + typedef HPolytope Hpolytope; + Hpolytope P; + + std::cout << "--- Testing volume of H-prod_simplex5" << std::endl; + P = gen_prod_simplex(5); + test_volume(P, std::pow(1.0 / factorial(5.0), 2)); + + std::cout << "--- Testing volume of H-prod_simplex10" << std::endl; + P = gen_prod_simplex(10); + test_volume(P, std::pow(1.0 / factorial(10.0), 2)); + + std::cout << "--- Testing volume of H-prod_simplex15" << std::endl; + P = gen_prod_simplex(15); + test_volume(P, std::pow(1.0 / factorial(15.0), 2)); + + std::cout << "--- Testing volume of H-prod_simplex20" << std::endl; + P = gen_prod_simplex(20); + test_volume(P, std::pow(1.0 / factorial(20.0), 2)); +} + +template +void call_test_simplex() { + typedef Cartesian Kernel; + typedef typename Kernel::Point Point; + typedef boost::mt19937 RNGType; + typedef HPolytope Hpolytope; + Hpolytope P; + + std::cout << "--- Testing volume of H-simplex10" << std::endl; + P = gen_simplex(10, false); + test_volume(P, 1.0 / factorial(10.0)); + + std::cout << "--- Testing volume of H-simplex20" << std::endl; + P = gen_simplex(20, false); + test_volume(P, 1.0 / factorial(20.0)); + + std::cout << "--- Testing volume of H-simplex30" << std::endl; + P = gen_simplex(30, false); + test_volume(P, 1.0 / factorial(30.0)); + + std::cout << "--- Testing volume of H-simplex40" << std::endl; + P = gen_simplex(40, false); + test_volume(P, 1.0 / factorial(40.0)); + + //std::cout << "--- Testing volume of H-simplex50" << std::endl; + //P = gen_simplex(50, false); + //test_volume(P, 1.0 / factorial(50.0)); +} + +template +void call_test_skinny_cube() { + typedef Cartesian Kernel; + typedef typename Kernel::Point Point; + typedef boost::mt19937 RNGType; + typedef HPolytope Hpolytope; + Hpolytope P; + + std::cout << "--- Testing volume of H-skinny_cube10" << std::endl; + P = gen_skinny_cube(10); + test_volume(P, 102400.0); + + //std::cout << "--- Testing volume of H-skinny_cube20" << std::endl; + //P = gen_skinny_cube(20); + //test_volume(P, 104857600.0); +} + + +TEST_CASE("cube") { + call_test_cube(); + //call_test_cube(); + //call_test_cube(); +} + +TEST_CASE("cross") { + call_test_cross(); + //call_test_cross(); + //call_test_cross(); +} + +TEST_CASE("birk") { + call_test_birk(); + //call_test_birk(); + //call_test_birk(); +} + +TEST_CASE("prod_simplex") { + call_test_prod_simplex(); + //call_test_prod_simplex(); + //call_test_prod_simplex(); +} + +TEST_CASE("simplex") { + call_test_simplex(); + //call_test_simplex(); + //call_test_simplex(); +} + +TEST_CASE("skinny_cube") { + call_test_skinny_cube(); + //call_test_skinny_cube(); + //call_test_skinny_cube(); +} From 43bc15a8a3d8bdb738044c92ca0b4cd70435bc56 Mon Sep 17 00:00:00 2001 From: Repouskos Panagiotis Date: Sat, 15 Feb 2020 17:05:20 +0200 Subject: [PATCH 08/21] Bug Error in creating Point --- include/convex_bodies/vpolytope.h | 2 +- include/samplers/samplers.h | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/include/convex_bodies/vpolytope.h b/include/convex_bodies/vpolytope.h index 0d2cb5393..96d3c8710 100644 --- a/include/convex_bodies/vpolytope.h +++ b/include/convex_bodies/vpolytope.h @@ -460,7 +460,7 @@ class VPolytope{ if (a.dot(V.row(outvert)) > 1.0) a = -a; a /= a.norm(); - Point s(_d, std::vector(&a[0], a.data()+a.cols()*a.rows())); + Point s(a); s = ((-2.0 * v.dot(s)) * s); v = s + v; diff --git a/include/samplers/samplers.h b/include/samplers/samplers.h index c479eb9ea..2b54affc4 100644 --- a/include/samplers/samplers.h +++ b/include/samplers/samplers.h @@ -483,7 +483,7 @@ void billiard_walk(ConvexBody &P, Point &p, NT diameter, VT &Ar, VT &Av, NT &lam return; } lambda_prev = dl * pbpair.first; - p = (lambda_prev * v) + p; + p += (lambda_prev * v); T -= lambda_prev; P.compute_reflection(v, p, pbpair.second); } @@ -498,7 +498,7 @@ void billiard_walk(ConvexBody &P, Point &p, NT diameter, VT &Ar, VT &Av, NT &lam } lambda_prev = dl * pbpair.first; - p = (lambda_prev * v) + p; + p += (lambda_prev * v); T -= lambda_prev; P.compute_reflection(v, p, pbpair.second); it++; From d589a1e307ba27445e7d3dad37b3fae2fdc114dd Mon Sep 17 00:00:00 2001 From: Repouskos Panagiotis Date: Sat, 15 Feb 2020 19:39:17 +0200 Subject: [PATCH 09/21] bug --- include/convex_bodies/zpolytope.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/convex_bodies/zpolytope.h b/include/convex_bodies/zpolytope.h index f074d58c3..5677d4406 100644 --- a/include/convex_bodies/zpolytope.h +++ b/include/convex_bodies/zpolytope.h @@ -389,7 +389,7 @@ class Zonotope { a = a/a.norm(); - Point s(_d, std::vector(&a[0], a.data()+a.cols()*a.rows())); + Point s(a); s = ((-2.0 * v.dot(s)) * s); v = s + v; From 755924000971acde07e10bab69605570847d13a3 Mon Sep 17 00:00:00 2001 From: Repouskos Panagiotis Date: Sat, 15 Feb 2020 19:39:59 +0200 Subject: [PATCH 10/21] use cooling balls in test --- test/volume_test_billiard.cpp | 37 ++++++++++++++++++++++++++++++++--- 1 file changed, 34 insertions(+), 3 deletions(-) diff --git a/test/volume_test_billiard.cpp b/test/volume_test_billiard.cpp index 1a4f7570c..e64972df2 100644 --- a/test/volume_test_billiard.cpp +++ b/test/volume_test_billiard.cpp @@ -5,18 +5,45 @@ // Licensed under GNU LGPL.3, see LICENCE file +#include "Eigen/Eigen" +#define VOLESTI_DEBUG +#include +#include "random.hpp" +#include "random/uniform_int.hpp" +#include "random/normal_distribution.hpp" +#include "random/uniform_real_distribution.hpp" +#include "volume.h" +#include "rotating.h" +#include "misc.h" +#include "linear_extensions.h" +#include "cooling_balls.h" +#include "cooling_hpoly.h" +#include "sample_only.h" +#include "exact_vols.h" +#include "Eigen/Eigen" +#include +#include "random.hpp" +#include "random/uniform_int.hpp" +#include "random/normal_distribution.hpp" +#include "random/uniform_real_distribution.hpp" +#include "volume.h" +#include "rotating.h" +#include "misc.h" #include "doctest.h" #include -#include "Eigen/Eigen" #include #include "random.hpp" #include "random/uniform_int.hpp" #include "random/normal_distribution.hpp" #include "random/uniform_real_distribution.hpp" #include "volume.h" +#include "cooling_hpoly.h" +#include "cooling_balls.h" #include "misc.h" #include "polytope_generators.h" #include +#include "sample_only.h" +#include "exact_vols.h" template NT factorial(NT n) @@ -32,7 +59,7 @@ void test_volume(Polytope &HP, NT expected, NT tolerance=0.1) // Setup the parameters int n = HP.dimension(); - int walk_len=10 + n/10; + int walk_len=1; int nexp=1, n_threads=1; NT e=1, err=0.0000000001; int rnum = std::pow(e,-2) * 400 * n * std::log(n); @@ -56,11 +83,15 @@ void test_volume(Polytope &HP, NT expected, NT tolerance=0.1) // Estimate the volume std::cout << "Number type: " << typeid(NT).name() << std::endl; NT vol = 0; + NT ball_radius=0.0, lb = 0.1, ub = 0.15, p = 0.75, rmax = 0.0, alpha = 0.2, round_val = 1.0; + NT C=2.0,ratio,frac=0.1,delta=-1.0; + vars_ban var_ban(lb, ub, p, rmax, alpha, 150, 125, 10, false); + unsigned int const num_of_exp = 10; for (unsigned int i=0; i Date: Fri, 21 Feb 2020 14:36:09 +0200 Subject: [PATCH 13/21] Optimizations - Use eigen function noalias() to avoid creating temporary copies when multiplying matrices - in getDirection() (samplers.h) don't use std::vector, only Eigen::Vector - fix bug in point.h, in constructor --- include/cartesian_geom/point.h | 9 ++++++--- include/convex_bodies/hpolytope.h | 23 +++++++++++------------ include/samplers/samplers.h | 16 +++++++++------- 3 files changed, 26 insertions(+), 22 deletions(-) diff --git a/include/cartesian_geom/point.h b/include/cartesian_geom/point.h index 3464e8b76..0f50b6c98 100644 --- a/include/cartesian_geom/point.h +++ b/include/cartesian_geom/point.h @@ -48,6 +48,7 @@ class point point(const unsigned int dim, std::vector cofs) { d = dim; + coeffs.resize(d); iter it = cofs.begin(); int i=0; @@ -80,6 +81,10 @@ class point return coeffs(i); } + FT* pointerToData() { + return coeffs.data(); + } + point operator+ (const point& p) const { point temp; temp.d = d; @@ -87,6 +92,7 @@ class point return temp; } + void operator+= (const point& p) { coeffs += p.getCoefficients(); } @@ -128,8 +134,6 @@ class point coeffs /= k; } - - bool operator== (point& p) const { int i=0; const Coeff & coeffs = p.getCoefficients(); @@ -143,7 +147,6 @@ class point return true; } - FT dot(const point& p) const { return coeffs.dot(p.getCoefficients()); } diff --git a/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h index d935a0505..d5705b66a 100644 --- a/include/convex_bodies/hpolytope.h +++ b/include/convex_bodies/hpolytope.h @@ -285,8 +285,8 @@ class HPolytope{ int m = num_of_hyperplanes(); - sum_nom = b - A * r.getCoefficients(); - sum_denom = A * v.getCoefficients(); + sum_nom.noalias() = b - A * r.getCoefficients(); + sum_denom.noalias() = A * v.getCoefficients(); NT* sum_nom_data = sum_nom.data(); NT* sum_denom_data = sum_denom.data(); @@ -321,9 +321,9 @@ class HPolytope{ unsigned int j; int m = num_of_hyperplanes(), facet; - Ar = A * r.getCoefficients(); + Ar.noalias() = A * r.getCoefficients(); sum_nom = b - Ar; - Av = A * v.getCoefficients();; + Av.noalias() = A * v.getCoefficients();; NT* Av_data = Av.data(); @@ -358,9 +358,9 @@ class HPolytope{ unsigned int j; int m = num_of_hyperplanes(), facet; - Ar += lambda_prev*Av; + Ar.noalias() += lambda_prev*Av; sum_nom = b - Ar; - Av = A * v.getCoefficients(); + Av.noalias() = A * v.getCoefficients(); NT* sum_nom_data = sum_nom.data(); NT* Av_data = Av.data(); @@ -411,7 +411,7 @@ class HPolytope{ viterator rit; sum_denom = A.col(rand_coord); - lamdas = b - A * r.getCoefficients(); + lamdas.noalias() = b - A * r.getCoefficients(); NT* lamda_data = lamdas.data(); NT* sum_denom_data = sum_denom.data(); @@ -445,7 +445,7 @@ class HPolytope{ int m = num_of_hyperplanes(); - lamdas += A.col(rand_coord_prev)* (r_prev[rand_coord_prev] - r[rand_coord_prev]); + lamdas.noalias() += A.col(rand_coord_prev)* (r_prev[rand_coord_prev] - r[rand_coord_prev]); NT* data = lamdas.data(); for (int i = 0; i < m; i++) { @@ -515,10 +515,9 @@ class HPolytope{ // VT a = A.row(facet); // Point s(_d, std::vector(&a[0], a.data()+a.cols()*a.rows())); - Point s(A.row(facet)); - s *= (-2.0 * v.dot(s)); - v += s; - +// Point s(A.row(facet)); +// s *= (-2.0 * v.dot(s)); + v += -2 * v.dot(A.row(facet)) * A.row(facet); } void free_them_all() {} diff --git a/include/samplers/samplers.h b/include/samplers/samplers.h index 2b54affc4..99e07391c 100644 --- a/include/samplers/samplers.h +++ b/include/samplers/samplers.h @@ -18,21 +18,23 @@ template Point get_direction(const unsigned int dim) { boost::normal_distribution<> rdist(0,1); - std::vector Xs(dim,0); NT normal = NT(0); unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); RNGType rng(seed); + + Point p(dim); + NT* data = p.pointerToData(); + //RNGType rng2 = var.rng; for (unsigned int i=0; i Date: Fri, 13 Mar 2020 20:07:01 +0200 Subject: [PATCH 14/21] leftovers from merge - make changes in code from last merge - delete tests I had previously added --- R-proj/src/sample_points.cpp | 12 +- include/generators/h_polytopes_gen.h | 7 +- include/volume/cooling_balls.h | 2 +- test/CMakeLists.txt | 4 - test/volume_test_billiard.cpp | 292 --------------------------- test/volume_test_rdhr.cpp | 247 ---------------------- 6 files changed, 8 insertions(+), 556 deletions(-) delete mode 100644 test/volume_test_billiard.cpp delete mode 100644 test/volume_test_rdhr.cpp diff --git a/R-proj/src/sample_points.cpp b/R-proj/src/sample_points.cpp index 0232ab38c..fe545601b 100644 --- a/R-proj/src/sample_points.cpp +++ b/R-proj/src/sample_points.cpp @@ -277,7 +277,7 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P = R_NilValue HP.normalize(); if (gaussian) { shift = MeanPoint; - HP.shift(Eigen::Map(&MeanPoint.get_coeffs()[0], MeanPoint.dimension())); + HP.shift(MeanPoint.getCoefficients()); MeanPoint = Point(dim); } break; @@ -296,7 +296,7 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P = R_NilValue if (billiard && diam < 0.0) VP.comp_diam(diam, 0.0); if (gaussian) { shift = MeanPoint; - VP.shift(Eigen::Map(&MeanPoint.get_coeffs()[0], MeanPoint.dimension())); + VP.shift(MeanPoint.getCoefficients()); MeanPoint = Point(dim); } break; @@ -315,7 +315,7 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P = R_NilValue if (billiard && diam < 0.0) ZP.comp_diam(diam, 0.0); if (gaussian) { shift = MeanPoint; - ZP.shift(Eigen::Map(&MeanPoint.get_coeffs()[0], MeanPoint.dimension())); + ZP.shift(MeanPoint.getCoefficients()); MeanPoint = Point(dim); } break; @@ -340,7 +340,7 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P = R_NilValue } if (gaussian) { shift = MeanPoint; - VPcVP.shift(Eigen::Map(&MeanPoint.get_coeffs()[0], MeanPoint.dimension())); + VPcVP.shift(MeanPoint.getCoefficients()); MeanPoint = Point(dim); } break; @@ -396,9 +396,9 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P = R_NilValue for (typename std::list::iterator rpit = randPoints.begin(); rpit!=randPoints.end(); rpit++, jj++) { if (gaussian) { - RetMat.col(jj) = Eigen::Map(&(*rpit).get_coeffs()[0], (*rpit).dimension()) + Eigen::Map(&shift.get_coeffs()[0], shift.dimension()); + RetMat.col(jj) = (*rpit).getCoefficients(); } else { - RetMat.col(jj) = Eigen::Map(&(*rpit).get_coeffs()[0], (*rpit).dimension()); + RetMat.col(jj) = (*rpit).getCoefficients(); } } diff --git a/include/generators/h_polytopes_gen.h b/include/generators/h_polytopes_gen.h index f0e273d64..b0319860a 100644 --- a/include/generators/h_polytopes_gen.h +++ b/include/generators/h_polytopes_gen.h @@ -34,13 +34,8 @@ Polytope random_hpoly(unsigned int dim, unsigned int m, double seed = std::numer for(unsigned int i=0; i(dim); - pit = p.iter_begin(); - j = 0; - for ( ; pit!=p.iter_end(); ++pit, ++j) { - A(i,j) = *pit; - } + A.row(i) = p.getCoefficients(); b(i) = 10.0; - } Polytope HP; HP.init(dim, A, b); diff --git a/include/volume/cooling_balls.h b/include/volume/cooling_balls.h index 06434a2af..56c7026c1 100644 --- a/include/volume/cooling_balls.h +++ b/include/volume/cooling_balls.h @@ -66,7 +66,7 @@ NT vol_cooling_balls(Polytope &P, UParameters &var, AParameters &var_ban, std::p var.che_rad = radius; // Move the chebychev center to the origin and apply the same shifting to the polytope - VT c_e = Eigen::Map(&c.get_coeffs()[0], c.dimension()); + VT c_e = c.getCoefficients(); P.shift(c_e); diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index d118fecb1..d93d0d185 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -73,8 +73,6 @@ else () add_library(test_main OBJECT test_main.cpp) add_executable (volume_test volume_test.cpp $) - add_executable (volume_test_billiard volume_test_billiard.cpp $) - add_executable (volume_test_rdhr volume_test_rdhr.cpp $) add_executable (cheb_test chebychev_test.cpp $) #add_executable (rounding_test rounding_test.cpp $) @@ -132,8 +130,6 @@ else () #TARGET_LINK_LIBRARIES(volume ${LP_SOLVE}) TARGET_LINK_LIBRARIES(generate ${LP_SOLVE}) TARGET_LINK_LIBRARIES(volume_test ${LP_SOLVE}) - TARGET_LINK_LIBRARIES(volume_test_rdhr ${LP_SOLVE}) - TARGET_LINK_LIBRARIES(volume_test_billiard ${LP_SOLVE}) TARGET_LINK_LIBRARIES(cheb_test ${LP_SOLVE}) #TARGET_LINK_LIBRARIES(rounding_test ${LP_SOLVE}) TARGET_LINK_LIBRARIES(volumeCG_test ${LP_SOLVE}) diff --git a/test/volume_test_billiard.cpp b/test/volume_test_billiard.cpp deleted file mode 100644 index a9e2ba8cc..000000000 --- a/test/volume_test_billiard.cpp +++ /dev/null @@ -1,292 +0,0 @@ -// VolEsti (volume computation and sampling library) - -// Copyright (c) 20012-2018 Vissarion Fisikopoulos -// Copyright (c) 2018 Apostolos Chalkis - -// Licensed under GNU LGPL.3, see LICENCE file - -#include "Eigen/Eigen" -#define VOLESTI_DEBUG -#include -#include "random.hpp" -#include "random/uniform_int.hpp" -#include "random/normal_distribution.hpp" -#include "random/uniform_real_distribution.hpp" -#include "volume.h" -#include "rotating.h" -#include "misc.h" -#include "linear_extensions.h" -#include "cooling_balls.h" -#include "cooling_hpoly.h" -#include "sample_only.h" -#include "exact_vols.h" -#include "Eigen/Eigen" -#include -#include "random.hpp" -#include "random/uniform_int.hpp" -#include "random/normal_distribution.hpp" -#include "random/uniform_real_distribution.hpp" -#include "volume.h" -#include "rotating.h" -#include "misc.h" -#include "doctest.h" -#include -#include -#include "random.hpp" -#include "random/uniform_int.hpp" -#include "random/normal_distribution.hpp" -#include "random/uniform_real_distribution.hpp" -#include "volume.h" -#include "cooling_hpoly.h" -#include "cooling_balls.h" -#include "misc.h" -#include "polytope_generators.h" -#include -#include "sample_only.h" -#include "exact_vols.h" - -template -NT factorial(NT n) -{ - return (n == 1 || n == 0) ? 1 : factorial(n - 1) * n; -} - -template -void test_volume(Polytope &HP, NT expected, NT tolerance=0.1) -{ - - typedef typename Polytope::PolytopePoint Point; - - // Setup the parameters - int n = HP.dimension(); - int walk_len=1; - int nexp=1, n_threads=1; - NT e=0.1, err=0.0000000001; - int rnum = std::pow(e,-2) * 400 * n * std::log(n); - unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); - RNGType rng(seed); - - boost::random::uniform_real_distribution<>(urdist); - boost::random::uniform_real_distribution<> urdist1(-1,1); - - - - std::pair InnerBall; - InnerBall = HP.ComputeInnerBall(); - NT diameter; - HP.comp_diam(diameter, InnerBall.second); - - - - //Compute chebychev ball// - std::pair CheBall; - - // Estimate the volume - std::cout << "Number type: " << typeid(NT).name() << std::endl; - NT vol = 0; - NT ball_radius=0.0, lb = 0.1, ub = 0.15, p = 0.75, rmax = 0.0, alpha = 0.2, round_val = 1.0; - NT C=2.0,ratio,frac=0.1,delta=-1.0; - - vars var(rnum,n,walk_len,n_threads,err,e,0,0.0,0,InnerBall.second,diameter,rng, - urdist,urdist1,-1.0,false,false,false,false,false,false,false, true,true); - vars_ban var_ban(lb, ub, p, rmax, alpha, 150, 125, 10, false); - - Polytope HP2 = HP; - unsigned int const num_of_exp = 10; - for (unsigned int i=0; i(10, false); - test_volume(P, 1024.0); - - std::cout << "--- Testing volume of H-cube20" << std::endl; - P = gen_cube(20, false); - test_volume(P, 1048576.0); - - std::cout << "--- Testing volume of H-cube30" << std::endl; - P = gen_cube(30, false); - test_volume(P, 1073742000.0, 0.2); -} - -template -void call_test_cross(){ - typedef Cartesian Kernel; - typedef typename Kernel::Point Point; - typedef boost::mt19937 RNGType; - typedef HPolytope Hpolytope; - - std::cout << "--- Testing volume of H-cross10" << std::endl; - Hpolytope P = gen_cross(10, false); - test_volume(P, 0.0002821869); -} - -template -void call_test_birk() { - typedef Cartesian Kernel; - typedef typename Kernel::Point Point; - typedef boost::mt19937 RNGType; - typedef HPolytope Hpolytope; - Hpolytope P; - - std::cout << "--- Testing volume of H-birk3" << std::endl; - std::ifstream inp; - std::vector > Pin; - inp.open("../R-proj/inst/extdata/birk3.ine",std::ifstream::in); - read_pointset(inp,Pin); - P.init(Pin); - test_volume(P, 0.125); - - std::cout << "--- Testing volume of H-birk4" << std::endl; - std::ifstream inp2; - std::vector > Pin2; - inp2.open("../R-proj/inst/extdata/birk4.ine",std::ifstream::in); - read_pointset(inp2,Pin2); - P.init(Pin2); - test_volume(P, 0.000970018); - - std::cout << "--- Testing volume of H-birk5" << std::endl; - std::ifstream inp3; - std::vector > Pin3; - inp3.open("../R-proj/inst/extdata/birk5.ine",std::ifstream::in); - read_pointset(inp3,Pin3); - P.init(Pin3); - test_volume(P, 0.000000225); - - std::cout << "--- Testing volume of H-birk6" << std::endl; - std::ifstream inp4; - std::vector > Pin4; - inp4.open("../R-proj/inst/extdata/birk6.ine",std::ifstream::in); - read_pointset(inp4,Pin4); - P.init(Pin4); - test_volume(P, 0.0000000000009455459196, 0.5); -} - -template -void call_test_prod_simplex() { - typedef Cartesian Kernel; - typedef typename Kernel::Point Point; - typedef boost::mt19937 RNGType; - typedef HPolytope Hpolytope; - Hpolytope P; - - std::cout << "--- Testing volume of H-prod_simplex5" << std::endl; - P = gen_prod_simplex(5); - test_volume(P, std::pow(1.0 / factorial(5.0), 2)); - - std::cout << "--- Testing volume of H-prod_simplex10" << std::endl; - P = gen_prod_simplex(10); - test_volume(P, std::pow(1.0 / factorial(10.0), 2)); - - std::cout << "--- Testing volume of H-prod_simplex15" << std::endl; - P = gen_prod_simplex(15); - test_volume(P, std::pow(1.0 / factorial(15.0), 2)); - - std::cout << "--- Testing volume of H-prod_simplex20" << std::endl; - P = gen_prod_simplex(20); - test_volume(P, std::pow(1.0 / factorial(20.0), 2)); -} - -template -void call_test_simplex() { - typedef Cartesian Kernel; - typedef typename Kernel::Point Point; - typedef boost::mt19937 RNGType; - typedef HPolytope Hpolytope; - Hpolytope P; - - std::cout << "--- Testing volume of H-simplex10" << std::endl; - P = gen_simplex(10, false); - test_volume(P, 1.0 / factorial(10.0)); - - std::cout << "--- Testing volume of H-simplex20" << std::endl; - P = gen_simplex(20, false); - test_volume(P, 1.0 / factorial(20.0)); - - std::cout << "--- Testing volume of H-simplex30" << std::endl; - P = gen_simplex(30, false); - test_volume(P, 1.0 / factorial(30.0)); - - std::cout << "--- Testing volume of H-simplex40" << std::endl; - P = gen_simplex(40, false); - test_volume(P, 1.0 / factorial(40.0)); - - //std::cout << "--- Testing volume of H-simplex50" << std::endl; - //P = gen_simplex(50, false); - //test_volume(P, 1.0 / factorial(50.0)); -} - -template -void call_test_skinny_cube() { - typedef Cartesian Kernel; - typedef typename Kernel::Point Point; - typedef boost::mt19937 RNGType; - typedef HPolytope Hpolytope; - Hpolytope P; - - std::cout << "--- Testing volume of H-skinny_cube10" << std::endl; - P = gen_skinny_cube(10); - test_volume(P, 102400.0); - - //std::cout << "--- Testing volume of H-skinny_cube20" << std::endl; - //P = gen_skinny_cube(20); - //test_volume(P, 104857600.0); -} - - -TEST_CASE("cube") { - call_test_cube(); - //call_test_cube(); - //call_test_cube(); -} - -TEST_CASE("cross") { - call_test_cross(); - //call_test_cross(); - //call_test_cross(); -} - -TEST_CASE("birk") { - call_test_birk(); - //call_test_birk(); - //call_test_birk(); -} - -TEST_CASE("prod_simplex") { - call_test_prod_simplex(); - //call_test_prod_simplex(); - //call_test_prod_simplex(); -} - -TEST_CASE("simplex") { - call_test_simplex(); - //call_test_simplex(); - //call_test_simplex(); -} - -TEST_CASE("skinny_cube") { - call_test_skinny_cube(); - //call_test_skinny_cube(); - //call_test_skinny_cube(); -} diff --git a/test/volume_test_rdhr.cpp b/test/volume_test_rdhr.cpp deleted file mode 100644 index 15ef868e3..000000000 --- a/test/volume_test_rdhr.cpp +++ /dev/null @@ -1,247 +0,0 @@ -// VolEsti (volume computation and sampling library) - -// Copyright (c) 20012-2018 Vissarion Fisikopoulos -// Copyright (c) 2018 Apostolos Chalkis - -// Licensed under GNU LGPL.3, see LICENCE file - -#include "doctest.h" -#include -#include "Eigen/Eigen" -#include -#include "random.hpp" -#include "random/uniform_int.hpp" -#include "random/normal_distribution.hpp" -#include "random/uniform_real_distribution.hpp" -#include "volume.h" -#include "misc.h" -#include "polytope_generators.h" -#include - -template -NT factorial(NT n) -{ - return (n == 1 || n == 0) ? 1 : factorial(n - 1) * n; -} - -template -void test_volume(Polytope &HP, NT expected, NT tolerance=0.1) -{ - - typedef typename Polytope::PolytopePoint Point; - - // Setup the parameters - int n = HP.dimension(); - int walk_len=10 + n/10; - int nexp=1, n_threads=1; - NT e=1, err=0.0000000001; - int rnum = std::pow(e,-2) * 400 * n * std::log(n); - unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); - RNGType rng(seed); - boost::normal_distribution<> rdist(0,1); - boost::random::uniform_real_distribution<>(urdist); - boost::random::uniform_real_distribution<> urdist1(-1,1); - - vars var(rnum,n,walk_len,n_threads,err,e,0,0,0,0,0.0,rng, - urdist,urdist1,-1.0,false,false,false,false,false,false,true, false,false); - - //Compute chebychev ball// - std::pair CheBall; - - // Estimate the volume - std::cout << "Number type: " << typeid(NT).name() << std::endl; - NT vol = 0; - unsigned int const num_of_exp = 10; - for (unsigned int i=0; i(10, false); - test_volume(P, 1024.0); - - std::cout << "--- Testing volume of H-cube20" << std::endl; - P = gen_cube(20, false); - test_volume(P, 1048576.0); - - std::cout << "--- Testing volume of H-cube30" << std::endl; - P = gen_cube(30, false); - test_volume(P, 1073742000.0, 0.2); -} - -template -void call_test_cross(){ - typedef Cartesian Kernel; - typedef typename Kernel::Point Point; - typedef boost::mt19937 RNGType; - typedef HPolytope Hpolytope; - - std::cout << "--- Testing volume of H-cross10" << std::endl; - Hpolytope P = gen_cross(10, false); - test_volume(P, 0.0002821869); -} - -template -void call_test_birk() { - typedef Cartesian Kernel; - typedef typename Kernel::Point Point; - typedef boost::mt19937 RNGType; - typedef HPolytope Hpolytope; - Hpolytope P; - - std::cout << "--- Testing volume of H-birk3" << std::endl; - std::ifstream inp; - std::vector > Pin; - inp.open("../R-proj/inst/extdata/birk3.ine",std::ifstream::in); - read_pointset(inp,Pin); - P.init(Pin); - test_volume(P, 0.125); - - std::cout << "--- Testing volume of H-birk4" << std::endl; - std::ifstream inp2; - std::vector > Pin2; - inp2.open("../R-proj/inst/extdata/birk4.ine",std::ifstream::in); - read_pointset(inp2,Pin2); - P.init(Pin2); - test_volume(P, 0.000970018); - - std::cout << "--- Testing volume of H-birk5" << std::endl; - std::ifstream inp3; - std::vector > Pin3; - inp3.open("../R-proj/inst/extdata/birk5.ine",std::ifstream::in); - read_pointset(inp3,Pin3); - P.init(Pin3); - test_volume(P, 0.000000225); - - std::cout << "--- Testing volume of H-birk6" << std::endl; - std::ifstream inp4; - std::vector > Pin4; - inp4.open("../R-proj/inst/extdata/birk6.ine",std::ifstream::in); - read_pointset(inp4,Pin4); - P.init(Pin4); - test_volume(P, 0.0000000000009455459196, 0.5); -} - -template -void call_test_prod_simplex() { - typedef Cartesian Kernel; - typedef typename Kernel::Point Point; - typedef boost::mt19937 RNGType; - typedef HPolytope Hpolytope; - Hpolytope P; - - std::cout << "--- Testing volume of H-prod_simplex5" << std::endl; - P = gen_prod_simplex(5); - test_volume(P, std::pow(1.0 / factorial(5.0), 2)); - - std::cout << "--- Testing volume of H-prod_simplex10" << std::endl; - P = gen_prod_simplex(10); - test_volume(P, std::pow(1.0 / factorial(10.0), 2)); - - std::cout << "--- Testing volume of H-prod_simplex15" << std::endl; - P = gen_prod_simplex(15); - test_volume(P, std::pow(1.0 / factorial(15.0), 2)); - - std::cout << "--- Testing volume of H-prod_simplex20" << std::endl; - P = gen_prod_simplex(20); - test_volume(P, std::pow(1.0 / factorial(20.0), 2)); -} - -template -void call_test_simplex() { - typedef Cartesian Kernel; - typedef typename Kernel::Point Point; - typedef boost::mt19937 RNGType; - typedef HPolytope Hpolytope; - Hpolytope P; - - std::cout << "--- Testing volume of H-simplex10" << std::endl; - P = gen_simplex(10, false); - test_volume(P, 1.0 / factorial(10.0)); - - std::cout << "--- Testing volume of H-simplex20" << std::endl; - P = gen_simplex(20, false); - test_volume(P, 1.0 / factorial(20.0)); - - std::cout << "--- Testing volume of H-simplex30" << std::endl; - P = gen_simplex(30, false); - test_volume(P, 1.0 / factorial(30.0)); - - std::cout << "--- Testing volume of H-simplex40" << std::endl; - P = gen_simplex(40, false); - test_volume(P, 1.0 / factorial(40.0)); - - //std::cout << "--- Testing volume of H-simplex50" << std::endl; - //P = gen_simplex(50, false); - //test_volume(P, 1.0 / factorial(50.0)); -} - -template -void call_test_skinny_cube() { - typedef Cartesian Kernel; - typedef typename Kernel::Point Point; - typedef boost::mt19937 RNGType; - typedef HPolytope Hpolytope; - Hpolytope P; - - std::cout << "--- Testing volume of H-skinny_cube10" << std::endl; - P = gen_skinny_cube(10); - test_volume(P, 102400.0); - - //std::cout << "--- Testing volume of H-skinny_cube20" << std::endl; - //P = gen_skinny_cube(20); - //test_volume(P, 104857600.0); -} - - -TEST_CASE("cube") { - call_test_cube(); - //call_test_cube(); - //call_test_cube(); -} - -TEST_CASE("cross") { - call_test_cross(); - //call_test_cross(); - //call_test_cross(); -} - -TEST_CASE("birk") { - call_test_birk(); - //call_test_birk(); - //call_test_birk(); -} - -TEST_CASE("prod_simplex") { - call_test_prod_simplex(); - //call_test_prod_simplex(); - //call_test_prod_simplex(); -} - -TEST_CASE("simplex") { - call_test_simplex(); - //call_test_simplex(); - //call_test_simplex(); -} - -TEST_CASE("skinny_cube") { - call_test_skinny_cube(); - //call_test_skinny_cube(); - //call_test_skinny_cube(); -} From a4e44a25927b60b23687b8a5d2caa841da1e4a5f Mon Sep 17 00:00:00 2001 From: Repouskos Panagiotis Date: Mon, 16 Mar 2020 22:11:35 +0200 Subject: [PATCH 15/21] cleanup - requested changes --- include/cartesian_geom/point.h | 5 ++--- include/convex_bodies/hpolytope.h | 8 ++++---- include/convex_bodies/zpolytope.h | 3 +-- include/rounding.h | 3 +-- include/samplers/samplers.h | 7 ++++--- 5 files changed, 12 insertions(+), 14 deletions(-) diff --git a/include/cartesian_geom/point.h b/include/cartesian_geom/point.h index 0f50b6c98..6c70873bf 100644 --- a/include/cartesian_geom/point.h +++ b/include/cartesian_geom/point.h @@ -18,11 +18,10 @@ class point private: unsigned int d; - typedef Eigen::Matrix Coeff; - - Coeff coeffs; + Eigen::Matrix coeffs; typedef typename std::vector::iterator iter; public: + typedef Eigen::Matrix Coeff; typedef typename K::FT FT; point() {} diff --git a/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h index 75fa7d884..57c221a6a 100644 --- a/include/convex_bodies/hpolytope.h +++ b/include/convex_bodies/hpolytope.h @@ -256,10 +256,11 @@ class HPolytope{ const NT* b_data = b.data(); for (int i = 0; i < m; i++) { - sum = *b_data - A.row(i) * p.getCoefficients(); - b_data++; //Check if corresponding hyperplane is violated - if (sum < NT(0)) return 0; + if (*b_data - A.row(i) * p.getCoefficients() < NT(0)) + return 0; + + b_data++; } return -1; } @@ -408,7 +409,6 @@ class HPolytope{ VT sum_denom; unsigned int j; int m = num_of_hyperplanes(); - viterator rit; sum_denom = A.col(rand_coord); lamdas.noalias() = b - A * r.getCoefficients(); diff --git a/include/convex_bodies/zpolytope.h b/include/convex_bodies/zpolytope.h index 135c7c9e8..49200aa5e 100644 --- a/include/convex_bodies/zpolytope.h +++ b/include/convex_bodies/zpolytope.h @@ -385,9 +385,8 @@ class Zonotope { } VT a = Fmat.fullPivLu().kernel(); - NT sum = p.getCoefficients().dot(a); - if(sum<0.0) a = -1.0*a; + if(p.getCoefficients().dot(a) < 0.0) a = -1.0*a; a = a/a.norm(); diff --git a/include/rounding.h b/include/rounding.h index 543d05011..3bea4ecff 100644 --- a/include/rounding.h +++ b/include/rounding.h @@ -16,7 +16,7 @@ #include "khach.h" -#include "../external/Eigen/Eigen" + // ----- ROUNDING ------ // // main rounding function @@ -26,7 +26,6 @@ std::pair rounding_min_ellipsoid(Polytope &P , const std::pair Coeff; unsigned int n=var.n, walk_len=var.walk_steps, i, j = 0; Point c = InnerBall.first; diff --git a/include/samplers/samplers.h b/include/samplers/samplers.h index 99e07391c..f0b50890e 100644 --- a/include/samplers/samplers.h +++ b/include/samplers/samplers.h @@ -115,7 +115,8 @@ void boundary_rand_point_generator(Polytope &P, { typedef typename Parameters::RNGType RNGType; typedef typename Point::FT NT; - typedef Eigen::Matrix VT; + typedef typename Point::Coeff VT; + unsigned int n = var.n; RNGType &rng = var.rng; boost::random::uniform_real_distribution<> urdist(0, 1); @@ -196,7 +197,7 @@ void rand_point_generator(Polytope &P, boost::random::uniform_real_distribution<> urdist(0, 1); boost::random::uniform_int_distribution<> uidist(0, n - 1); - typedef Eigen::Matrix VT; + typedef typename Point::Coeff VT; VT lamdas, Av; lamdas.setZero(P.num_of_hyperplanes()); @@ -267,7 +268,7 @@ void rand_point_generator(BallPoly &PBLarge, boost::random::uniform_real_distribution<> urdist(0, 1); boost::random::uniform_int_distribution<> uidist(0, n - 1); - typedef Eigen::Matrix VT; + typedef typename Point::Coeff VT; VT lamdas, Av; lamdas.setZero(PBLarge.num_of_hyperplanes()); From 623521e1b47f13d7a99429c82321e587ff9c6069 Mon Sep 17 00:00:00 2001 From: Repouskos Panagiotis Date: Mon, 16 Mar 2020 22:19:07 +0200 Subject: [PATCH 16/21] update copyrights --- include/cartesian_geom/point.h | 1 + include/convex_bodies/ball.h | 2 ++ include/convex_bodies/ballintersectconvex.h | 1 + include/convex_bodies/ellipsoids.h | 1 + include/convex_bodies/hpolytope.h | 1 + include/convex_bodies/vpolyintersectvpoly.h | 1 + include/convex_bodies/vpolytope.h | 1 + include/convex_bodies/zonoIntersecthpoly.h | 2 ++ include/convex_bodies/zpolytope.h | 1 + include/lp_oracles/misc_lp.h | 2 ++ include/lp_oracles/solve_lp.h | 1 + include/lp_oracles/vpolyoracles.h | 2 ++ include/lp_oracles/zpolyoracles.h | 2 ++ include/rotating.h | 2 ++ include/rounding.h | 1 + include/samplers/gaussian_samplers.h | 1 + include/samplers/sample_only.h | 1 + include/samplers/samplers.h | 1 + include/samplers/simplex_samplers.h | 1 + include/volume/cooling_balls.h | 2 ++ include/volume/cooling_hpoly.h | 2 ++ include/volume/copulas.h | 1 + include/volume/exact_vols.h | 2 ++ include/volume/volume.h | 1 + 24 files changed, 33 insertions(+) diff --git a/include/cartesian_geom/point.h b/include/cartesian_geom/point.h index 6c70873bf..b3708f1b3 100644 --- a/include/cartesian_geom/point.h +++ b/include/cartesian_geom/point.h @@ -3,6 +3,7 @@ // Copyright (c) 2018 Vissarion Fisikopoulos, Apostolos Chalkis //Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 program. +//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. // Licensed under GNU LGPL.3, see LICENCE file diff --git a/include/convex_bodies/ball.h b/include/convex_bodies/ball.h index e5c3b875e..0d2ceb36d 100644 --- a/include/convex_bodies/ball.h +++ b/include/convex_bodies/ball.h @@ -3,6 +3,8 @@ // Copyright (c) 20012-2020 Vissarion Fisikopoulos // Copyright (c) 2018-2020 Apostolos Chalkis +//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. + // Licensed under GNU LGPL.3, see LICENCE file #ifndef BALL_H diff --git a/include/convex_bodies/ballintersectconvex.h b/include/convex_bodies/ballintersectconvex.h index 288b84fb9..7cbe3e7e7 100644 --- a/include/convex_bodies/ballintersectconvex.h +++ b/include/convex_bodies/ballintersectconvex.h @@ -4,6 +4,7 @@ // Copyright (c) 2018-2020 Apostolos Chalkis //Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 program. +//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. // Licensed under GNU LGPL.3, see LICENCE file diff --git a/include/convex_bodies/ellipsoids.h b/include/convex_bodies/ellipsoids.h index 19eae94f8..72ef9cc88 100644 --- a/include/convex_bodies/ellipsoids.h +++ b/include/convex_bodies/ellipsoids.h @@ -4,6 +4,7 @@ // Copyright (c) 2018 Apostolos Chalkis //Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 program. +//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. // Licensed under GNU LGPL.3, see LICENCE file diff --git a/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h index 57c221a6a..bbc5162b9 100644 --- a/include/convex_bodies/hpolytope.h +++ b/include/convex_bodies/hpolytope.h @@ -4,6 +4,7 @@ // Copyright (c) 2018 Apostolos Chalkis //Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 program. +//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. // Licensed under GNU LGPL.3, see LICENCE file diff --git a/include/convex_bodies/vpolyintersectvpoly.h b/include/convex_bodies/vpolyintersectvpoly.h index 4b5b09a97..a74f6e863 100644 --- a/include/convex_bodies/vpolyintersectvpoly.h +++ b/include/convex_bodies/vpolyintersectvpoly.h @@ -4,6 +4,7 @@ // Copyright (c) 2018 Apostolos Chalkis //Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 program. +//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. // Licensed under GNU LGPL.3, see LICENCE file diff --git a/include/convex_bodies/vpolytope.h b/include/convex_bodies/vpolytope.h index 9207b5668..c841cd9f1 100644 --- a/include/convex_bodies/vpolytope.h +++ b/include/convex_bodies/vpolytope.h @@ -4,6 +4,7 @@ // Copyright (c) 2018 Apostolos Chalkis //Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 program. +//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. // Licensed under GNU LGPL.3, see LICENCE file diff --git a/include/convex_bodies/zonoIntersecthpoly.h b/include/convex_bodies/zonoIntersecthpoly.h index 08ead19d1..8ebae3148 100644 --- a/include/convex_bodies/zonoIntersecthpoly.h +++ b/include/convex_bodies/zonoIntersecthpoly.h @@ -3,6 +3,8 @@ // Copyright (c) 20012-2018 Vissarion Fisikopoulos // Copyright (c) 2018 Apostolos Chalkis +//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. + #ifndef ZONOINTERSECTHPOLY_H #define ZONOINTERSECTHPOLY_H diff --git a/include/convex_bodies/zpolytope.h b/include/convex_bodies/zpolytope.h index 49200aa5e..818a81759 100644 --- a/include/convex_bodies/zpolytope.h +++ b/include/convex_bodies/zpolytope.h @@ -4,6 +4,7 @@ // Copyright (c) 2018 Apostolos Chalkis //Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 program. +//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. // Licensed under GNU LGPL.3, see LICENCE file diff --git a/include/lp_oracles/misc_lp.h b/include/lp_oracles/misc_lp.h index c87baa174..fb78bf2aa 100644 --- a/include/lp_oracles/misc_lp.h +++ b/include/lp_oracles/misc_lp.h @@ -1,6 +1,8 @@ #ifndef MISC_LP_H #define MISC_LP_H +//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. + #include #include #include diff --git a/include/lp_oracles/solve_lp.h b/include/lp_oracles/solve_lp.h index a3224e67f..a5a079afb 100644 --- a/include/lp_oracles/solve_lp.h +++ b/include/lp_oracles/solve_lp.h @@ -3,6 +3,7 @@ // Copyright (c) 2018 Vissarion Fisikopoulos, Apostolos Chalkis //Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 program. +//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. // VolEsti is free software: you can redistribute it and/or modify it // under the terms of the GNU Lesser General Public License as published by diff --git a/include/lp_oracles/vpolyoracles.h b/include/lp_oracles/vpolyoracles.h index 9ae2e68fd..94c872928 100644 --- a/include/lp_oracles/vpolyoracles.h +++ b/include/lp_oracles/vpolyoracles.h @@ -1,6 +1,8 @@ #ifndef VPOLYORACLES_H #define VPOLYORACLES_H +//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. + #include #include #include diff --git a/include/lp_oracles/zpolyoracles.h b/include/lp_oracles/zpolyoracles.h index 3b8e853ef..1396e3688 100644 --- a/include/lp_oracles/zpolyoracles.h +++ b/include/lp_oracles/zpolyoracles.h @@ -1,6 +1,8 @@ #ifndef ZPOLYORACLES_H #define ZPOLYORACLES_H +//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. + #include #include #include diff --git a/include/rotating.h b/include/rotating.h index a29c705cb..c9e4efd25 100644 --- a/include/rotating.h +++ b/include/rotating.h @@ -3,6 +3,8 @@ // Copyright (c) 20012-2019 Vissarion Fisikopoulos // Copyright (c) 2019 Apostolos Chalkis +//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. + // Licensed under GNU LGPL.3, see LICENCE file #ifndef ROTATING_H diff --git a/include/rounding.h b/include/rounding.h index 3bea4ecff..f4c35e900 100644 --- a/include/rounding.h +++ b/include/rounding.h @@ -4,6 +4,7 @@ // Copyright (c) 2018 Apostolos Chalkis //Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 program. +//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. // The functions in this header file call Bojan Nikolic implementation // of Todd and Yildirim algorithm in "On Khachiyan's Algorithm for the Computation of Minimum Volume Enclosing Ellipsoids", 2005 diff --git a/include/samplers/gaussian_samplers.h b/include/samplers/gaussian_samplers.h index 0fbfdecf5..1982dbb01 100644 --- a/include/samplers/gaussian_samplers.h +++ b/include/samplers/gaussian_samplers.h @@ -3,6 +3,7 @@ // Copyright (c) 2018 Vissarion Fisikopoulos, Apostolos Chalkis //Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 program. +//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. // VolEsti is free software: you can redistribute it and/or modify it // under the terms of the GNU Lesser General Public License as published by diff --git a/include/samplers/sample_only.h b/include/samplers/sample_only.h index 7bbcacd00..2c241890d 100644 --- a/include/samplers/sample_only.h +++ b/include/samplers/sample_only.h @@ -3,6 +3,7 @@ // Copyright (c) 2018 Vissarion Fisikopoulos, Apostolos Chalkis //Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 program. +//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. // VolEsti is free software: you can redistribute it and/or modify it // under the terms of the GNU Lesser General Public License as published by diff --git a/include/samplers/samplers.h b/include/samplers/samplers.h index f0b50890e..c901b24ec 100644 --- a/include/samplers/samplers.h +++ b/include/samplers/samplers.h @@ -4,6 +4,7 @@ // Copyright (c) 2018 Apostolos Chalkis //Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 program. +//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. // Licensed under GNU LGPL.3, see LICENCE file diff --git a/include/samplers/simplex_samplers.h b/include/samplers/simplex_samplers.h index 28774d15c..563bdd410 100644 --- a/include/samplers/simplex_samplers.h +++ b/include/samplers/simplex_samplers.h @@ -3,6 +3,7 @@ // Copyright (c) 2018 Vissarion Fisikopoulos, Apostolos Chalkis //Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 program. +//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. // VolEsti is free software: you can redistribute it and/or modify it // under the terms of the GNU Lesser General Public License as published by diff --git a/include/volume/cooling_balls.h b/include/volume/cooling_balls.h index 56c7026c1..0adf461fa 100644 --- a/include/volume/cooling_balls.h +++ b/include/volume/cooling_balls.h @@ -3,6 +3,8 @@ // Copyright (c) 20012-2019 Vissarion Fisikopoulos // Copyright (c) 2018-2019 Apostolos Chalkis +//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. + #ifndef COOLING_BALLS_H #define COOLING_BALLS_H diff --git a/include/volume/cooling_hpoly.h b/include/volume/cooling_hpoly.h index 57777bd5a..f134d686d 100644 --- a/include/volume/cooling_hpoly.h +++ b/include/volume/cooling_hpoly.h @@ -3,6 +3,8 @@ // Copyright (c) 20012-2018 Vissarion Fisikopoulos // Copyright (c) 2018 Apostolos Chalkis +//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. + #ifndef COOLING_HPOLY_H #define COOLING_HPOLY_H diff --git a/include/volume/copulas.h b/include/volume/copulas.h index cf163ac54..0d2ec58b2 100644 --- a/include/volume/copulas.h +++ b/include/volume/copulas.h @@ -3,6 +3,7 @@ // Copyright (c) 2018 Vissarion Fisikopoulos, Apostolos Chalkis //Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 program. +//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. // VolEsti is free software: you can redistribute it and/or modify it // under the terms of the GNU Lesser General Public License as published by diff --git a/include/volume/exact_vols.h b/include/volume/exact_vols.h index 71525f63d..ceb97e18b 100644 --- a/include/volume/exact_vols.h +++ b/include/volume/exact_vols.h @@ -3,6 +3,8 @@ // Copyright (c) 20012-2018 Vissarion Fisikopoulos // Copyright (c) 2018 Apostolos Chalkis +//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. + // Licensed under GNU LGPL.3, see LICENCE file #ifndef ZONOTOPE_EXACT_VOL_H diff --git a/include/volume/volume.h b/include/volume/volume.h index 30d9ed434..72c49aa79 100644 --- a/include/volume/volume.h +++ b/include/volume/volume.h @@ -5,6 +5,7 @@ // Copyright (c) 2018 Vissarion Fisikopoulos, Apostolos Chalkis //Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 program. +//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. // Licensed under GNU LGPL.3, see LICENCE file From f3fbec0af5dd8ab8663537b9663fad72b2922b21 Mon Sep 17 00:00:00 2001 From: Repouskos Panagiotis Date: Mon, 16 Mar 2020 22:27:24 +0200 Subject: [PATCH 17/21] bug --- include/samplers/samplers.h | 1 - 1 file changed, 1 deletion(-) diff --git a/include/samplers/samplers.h b/include/samplers/samplers.h index c901b24ec..4dcf7994a 100644 --- a/include/samplers/samplers.h +++ b/include/samplers/samplers.h @@ -123,7 +123,6 @@ void boundary_rand_point_generator(Polytope &P, boost::random::uniform_real_distribution<> urdist(0, 1); boost::random::uniform_int_distribution<> uidist(0, n - 1); - typedef Eigen::Matrix VT; VT lamdas, Av; lamdas.setZero(P.num_of_hyperplanes()); From 7b2794a7bbed475edd7c23320ae17b98a80f306f Mon Sep 17 00:00:00 2001 From: Repouskos Panagiotis Date: Mon, 16 Mar 2020 23:19:55 +0200 Subject: [PATCH 18/21] bug --- include/cartesian_geom/point.h | 1 - 1 file changed, 1 deletion(-) diff --git a/include/cartesian_geom/point.h b/include/cartesian_geom/point.h index b3708f1b3..f27a67957 100644 --- a/include/cartesian_geom/point.h +++ b/include/cartesian_geom/point.h @@ -11,7 +11,6 @@ #define POINT_H #include -#include "../external/Eigen/Eigen" template class point From fc21da5274fe40eda55e1611bace21aaf89d9554 Mon Sep 17 00:00:00 2001 From: Repouskos Panagiotis Date: Sat, 21 Mar 2020 15:01:34 +0200 Subject: [PATCH 19/21] requested changes --- include/annealing/ratio_estimation.h | 1 + include/cartesian_geom/point.h | 16 +++++----------- include/convex_bodies/ball.h | 6 ++---- include/convex_bodies/hpolytope.h | 2 +- include/convex_bodies/vpolytope.h | 9 +++------ include/convex_bodies/zpolytope.h | 7 +++---- include/lp_oracles/misc_lp.h | 2 -- include/lp_oracles/vpolyoracles.h | 1 - include/lp_oracles/zpolyoracles.h | 1 - include/rotating.h | 1 - include/rounding.h | 10 ++++------ include/samplers/gaussian_samplers.h | 2 +- include/samplers/sample_only.h | 1 - include/samplers/samplers.h | 4 +--- include/samplers/simplex_samplers.h | 1 - include/volume/cooling_balls.h | 3 +-- include/volume/cooling_hpoly.h | 1 - include/volume/copulas.h | 1 - 18 files changed, 22 insertions(+), 47 deletions(-) diff --git a/include/annealing/ratio_estimation.h b/include/annealing/ratio_estimation.h index cb66508e3..d8256fb60 100644 --- a/include/annealing/ratio_estimation.h +++ b/include/annealing/ratio_estimation.h @@ -3,6 +3,7 @@ // Copyright (c) 20012-2018 Vissarion Fisikopoulos // Copyright (c) 2018 Apostolos Chalkis +//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. #ifndef RATIO_ESTIMATION_H #define RATIO_ESTIMATION_H diff --git a/include/cartesian_geom/point.h b/include/cartesian_geom/point.h index f27a67957..08f8b6593 100644 --- a/include/cartesian_geom/point.h +++ b/include/cartesian_geom/point.h @@ -12,7 +12,7 @@ #include -template +template class point { private: @@ -136,6 +136,9 @@ class point bool operator== (point& p) const { int i=0; const Coeff & coeffs = p.getCoefficients(); + + /* degree of approximation in + "The art of computer programming" (vol II), p. 234, Donald. E. Knuth. */ FT e = 0.00000000001; for (i=0 ; icoeffs(i) - coeffs(i)) > e *std::abs(this->coeffs(i)) || @@ -176,18 +179,9 @@ class point #endif } - -// iter iter_begin() { -// return coeffs.begin(); -// } -// -// iter iter_end() { -// return coeffs.end(); -// } - }; -template +template point operator* (const typename K::FT& k, point& p) { return p * k; } diff --git a/include/convex_bodies/ball.h b/include/convex_bodies/ball.h index 0d2ceb36d..91a7cfb82 100644 --- a/include/convex_bodies/ball.h +++ b/include/convex_bodies/ball.h @@ -57,14 +57,12 @@ struct Ball{ return std::pair ((NT(-1)*vrc + disc_sqrt)/v2, (NT(-1)*vrc - disc_sqrt)/v2); } - std::pair line_intersect(Point &r, Point &v, const VT &Ar, - const VT &Av){ + std::pair line_intersect(Point &r, Point &v, const VT &Ar, const VT &Av){ return line_intersect(r, v); } - std::pair line_intersect(Point &r, Point &v, const VT &Ar, - const VT &Av, NT &lambda_prev) { + std::pair line_intersect(Point &r, Point &v, const VT &Ar, const VT &Av, NT &lambda_prev) { return line_intersect(r, v); } diff --git a/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h index 4ed9d12d5..afb03b875 100644 --- a/include/convex_bodies/hpolytope.h +++ b/include/convex_bodies/hpolytope.h @@ -474,7 +474,7 @@ class HPolytope{ // shift polytope by a point c void shift(const VT &c){ - b = b - A*c; + b -= A*c; } diff --git a/include/convex_bodies/vpolytope.h b/include/convex_bodies/vpolytope.h index aa3e54291..14d0135f8 100644 --- a/include/convex_bodies/vpolytope.h +++ b/include/convex_bodies/vpolytope.h @@ -301,8 +301,6 @@ class VPolytope{ boost::numeric::ublas::matrix Ap(_d,randPoints.size()); typename std::list::iterator rpit=randPoints.begin(); - - unsigned int i, j = 0; for ( ; rpit!=randPoints.end(); rpit++, j++) { const NT* point_data = rpit->getCoefficients().data(); @@ -466,10 +464,9 @@ class VPolytope{ if (a.dot(V.row(outvert)) > 1.0) a = -a; a /= a.norm(); - Point s(a); - - s = ((-2.0 * v.dot(s)) * s); - v = s + v; + // compute reflection + a *= (-2.0 * v.dot(a)); + v += a; } void free_them_all() { diff --git a/include/convex_bodies/zpolytope.h b/include/convex_bodies/zpolytope.h index 5ee5b84ec..e79f5a122 100644 --- a/include/convex_bodies/zpolytope.h +++ b/include/convex_bodies/zpolytope.h @@ -392,10 +392,9 @@ class Zonotope { a = a/a.norm(); - Point s(a); - - s = ((-2.0 * v.dot(s)) * s); - v = s + v; + // compute reflection + a *= (-2.0 * v.dot(a)); + v += a; } void free_them_all() { diff --git a/include/lp_oracles/misc_lp.h b/include/lp_oracles/misc_lp.h index fb78bf2aa..c87baa174 100644 --- a/include/lp_oracles/misc_lp.h +++ b/include/lp_oracles/misc_lp.h @@ -1,8 +1,6 @@ #ifndef MISC_LP_H #define MISC_LP_H -//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. - #include #include #include diff --git a/include/lp_oracles/vpolyoracles.h b/include/lp_oracles/vpolyoracles.h index 94c872928..e988b055c 100644 --- a/include/lp_oracles/vpolyoracles.h +++ b/include/lp_oracles/vpolyoracles.h @@ -1,7 +1,6 @@ #ifndef VPOLYORACLES_H #define VPOLYORACLES_H -//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. #include #include diff --git a/include/lp_oracles/zpolyoracles.h b/include/lp_oracles/zpolyoracles.h index 1396e3688..06ad2c686 100644 --- a/include/lp_oracles/zpolyoracles.h +++ b/include/lp_oracles/zpolyoracles.h @@ -1,7 +1,6 @@ #ifndef ZPOLYORACLES_H #define ZPOLYORACLES_H -//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. #include #include diff --git a/include/rotating.h b/include/rotating.h index c9e4efd25..366b8c5f1 100644 --- a/include/rotating.h +++ b/include/rotating.h @@ -3,7 +3,6 @@ // Copyright (c) 20012-2019 Vissarion Fisikopoulos // Copyright (c) 2019 Apostolos Chalkis -//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. // Licensed under GNU LGPL.3, see LICENCE file diff --git a/include/rounding.h b/include/rounding.h index f4c35e900..63cd369f6 100644 --- a/include/rounding.h +++ b/include/rounding.h @@ -55,13 +55,12 @@ std::pair rounding_min_ellipsoid(Polytope &P , const std::pair Ap(n,randPoints.size()); typename std::list::iterator rpit=randPoints.begin(); - for ( ; rpit!=randPoints.end(); rpit++, j++) { for (i=0 ; idimension(); i++){ Ap(i,j)=double((*rpit)[i]); } } - boost::numeric::ublas::matrix Q(n,n); + boost::numeric::ublas::matrix Q(n,n); //TODO: remove dependence on ublas and copy to eigen boost::numeric::ublas::vector c2(n); size_t w=1000; KhachiyanAlgo(Ap,0.01,w,Q,c2); // call Khachiyan algorithm @@ -108,16 +107,15 @@ void get_vpoly_center(Polytope &P) { typedef typename Polytope::VT VT; typedef typename Polytope::PolytopePoint Point; - unsigned int n = P.dimension(), i, j = 0; + unsigned int n = P.dimension(); std::list randPoints; //ds for storing rand points P.get_points_for_rounding(randPoints); boost::numeric::ublas::matrix Ap(n,randPoints.size()); typename std::list::iterator rpit=randPoints.begin(); - for ( ; rpit!=randPoints.end(); rpit++, j++) { - i=0; - for ( ; idimension(); i++){ + for (int j=0 ; rpit!=randPoints.end(); rpit++, j++) { + for (int i=0 ; idimension(); i++){ Ap(i,j)=double((*rpit)[i]); } } diff --git a/include/samplers/gaussian_samplers.h b/include/samplers/gaussian_samplers.h index 1982dbb01..bde4ddd83 100644 --- a/include/samplers/gaussian_samplers.h +++ b/include/samplers/gaussian_samplers.h @@ -286,7 +286,7 @@ void gaussian_hit_and_run(Point &p, // hit-and-run with orthogonal directions and update -template +template void gaussian_hit_and_run_coord_update(Point &p, Point &p_prev, Polytope &P, diff --git a/include/samplers/sample_only.h b/include/samplers/sample_only.h index 454d29891..23034fc1d 100644 --- a/include/samplers/sample_only.h +++ b/include/samplers/sample_only.h @@ -3,7 +3,6 @@ // Copyright (c) 2018 Vissarion Fisikopoulos, Apostolos Chalkis //Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 program. -//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. // VolEsti is free software: you can redistribute it and/or modify it // under the terms of the GNU Lesser General Public License as published by diff --git a/include/samplers/samplers.h b/include/samplers/samplers.h index 4117f26a1..9754af0fb 100644 --- a/include/samplers/samplers.h +++ b/include/samplers/samplers.h @@ -27,7 +27,7 @@ Point get_direction(const unsigned int dim) { NT* data = p.pointerToData(); //RNGType rng2 = var.rng; - for (unsigned int i=0; i(P, BallSet, ratios, N * nu, nu, lb, ub, radius, alpha, var, rmax) ){ diff --git a/include/volume/cooling_hpoly.h b/include/volume/cooling_hpoly.h index f134d686d..7b0025900 100644 --- a/include/volume/cooling_hpoly.h +++ b/include/volume/cooling_hpoly.h @@ -3,7 +3,6 @@ // Copyright (c) 20012-2018 Vissarion Fisikopoulos // Copyright (c) 2018 Apostolos Chalkis -//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. #ifndef COOLING_HPOLY_H #define COOLING_HPOLY_H diff --git a/include/volume/copulas.h b/include/volume/copulas.h index 0d2ec58b2..cf163ac54 100644 --- a/include/volume/copulas.h +++ b/include/volume/copulas.h @@ -3,7 +3,6 @@ // Copyright (c) 2018 Vissarion Fisikopoulos, Apostolos Chalkis //Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 program. -//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. // VolEsti is free software: you can redistribute it and/or modify it // under the terms of the GNU Lesser General Public License as published by From 89655c3963fcc6b6c599d3cb2d2b2bbe66fec4d0 Mon Sep 17 00:00:00 2001 From: Repouskos Panagiotis Date: Mon, 23 Mar 2020 13:09:00 +0200 Subject: [PATCH 20/21] use += *= operators with points --- include/samplers/samplers.h | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/include/samplers/samplers.h b/include/samplers/samplers.h index 9754af0fb..97e40f186 100644 --- a/include/samplers/samplers.h +++ b/include/samplers/samplers.h @@ -60,7 +60,7 @@ Point get_point_in_Dsphere(const unsigned int dim, const NT &radius){ Point p = get_direction(dim); U = urdist(rng2); U = std::pow(U, 1.0/(NT(dim))); - p = (radius*U)*p; + p *= (radius*U); return p; } @@ -72,7 +72,7 @@ void ball_walk(Point &p, { //typedef typename Parameters::RNGType RNGType; Point y = get_point_in_Dsphere(p.dimension(), delta); - y = y + p; + y += p; if (P.is_in(y)==-1) p = y; } @@ -143,7 +143,7 @@ void boundary_rand_point_generator(Polytope &P, v = get_direction(n); std::pair bpair = P.line_intersect(p, v, lamdas, Av); lambda = urdist(rng) * (bpair.first - bpair.second) + bpair.second; - p = (lambda * v) + p; + p += (lambda * v); } //hit_and_run(p, P, var); @@ -169,7 +169,7 @@ void boundary_rand_point_generator(Polytope &P, p1 = (bpair.first * v) + p; p2 = (bpair.second * v) + p; lambda = urdist(rng) * (bpair.first - bpair.second) + bpair.second; - p = (lambda * v) + p; + p += (lambda * v); } } @@ -217,7 +217,7 @@ void rand_point_generator(Polytope &P, v = get_direction(n); std::pair bpair = P.line_intersect(p, v, lamdas, Av); lambda = urdist(rng) * (bpair.first - bpair.second) + bpair.second; - p = (lambda * v) + p; + p += (lambda * v); //hit_and_run(p, P, var); } else { billiard_walk(P, p, var.diameter, lamdas, Av, lambda, var, true); @@ -237,7 +237,7 @@ void rand_point_generator(Polytope &P, v = get_direction(n); std::pair bpair = P.line_intersect(p, v, lamdas, Av, lambda); lambda = urdist(rng) * (bpair.first - bpair.second) + bpair.second; - p = (lambda * v) + p; + p += (lambda * v); //hit_and_run(p, P, var); } else { billiard_walk(P, p, var.diameter, lamdas, Av, lambda, var); @@ -290,7 +290,7 @@ void rand_point_generator(BallPoly &PBLarge, v = get_direction(n); std::pair bpair = PBLarge.line_intersect(p, v, lamdas, Av); lambda = urdist(rng) * (bpair.first - bpair.second) + bpair.second; - p = (lambda * v) + p; + p += (lambda * v); //hit_and_run(p, PBLarge, var); } else { billiard_walk(PBLarge, p, var.diameter, lamdas, Av, lambda, var, true); @@ -309,7 +309,7 @@ void rand_point_generator(BallPoly &PBLarge, v = get_direction(n); std::pair bpair = PBLarge.line_intersect(p, v, lamdas, Av, lambda); lambda = urdist(rng) * (bpair.first - bpair.second) + bpair.second; - p = (lambda * v) + p; + p += (lambda * v); //hit_and_run(p, PBLarge, var); } else { billiard_walk(PBLarge, p, var.diameter, lamdas, Av, lambda, var); @@ -354,7 +354,7 @@ void uniform_first_point(Polytope &P, v = get_direction(n); std::pair bpair = P.line_intersect(p, v, lamdas, Av); lambda = urdist(rng) * (bpair.first - bpair.second) + bpair.second; - p = (lambda * v) + p; + p += (lambda * v); } else { billiard_walk(P, p, var.diameter, lamdas, Av, lambda, var, true); } @@ -377,7 +377,7 @@ void uniform_first_point(Polytope &P, v = get_direction(n); std::pair bpair = P.line_intersect(p, v, lamdas, Av, lambda); lambda = urdist(rng) * (bpair.first - bpair.second) + bpair.second; - p = (lambda * v) + p; + p += (lambda * v); } } else { billiard_walk(P, p, var.diameter, lamdas, Av, lambda, var); @@ -413,7 +413,7 @@ void uniform_next_point(Polytope &P, v = get_direction(n); std::pair bpair = P.line_intersect(p, v, lamdas, Av, lambda); lambda = urdist(rng) * (bpair.first - bpair.second) + bpair.second; - p = (lambda * v) + p; + p += (lambda * v); } } else if (var.bill_walk) { for (unsigned int j = 0; j < walk_len; j++) billiard_walk(P, p, var.diameter, lamdas, Av, lambda, var); From 09b82499c50f8fecd733b7721f298ba4e534f328 Mon Sep 17 00:00:00 2001 From: Repouskos Panagiotis Date: Mon, 23 Mar 2020 14:27:12 +0200 Subject: [PATCH 21/21] use += *= operators with points --- include/samplers/gaussian_samplers.h | 4 ++-- include/samplers/samplers.h | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/include/samplers/gaussian_samplers.h b/include/samplers/gaussian_samplers.h index bde4ddd83..607e240fe 100644 --- a/include/samplers/gaussian_samplers.h +++ b/include/samplers/gaussian_samplers.h @@ -60,7 +60,7 @@ NT get_max_coord(const NT &l, const NT &u, const NT &a_i) { // Pick a point from the distribution exp(-a_i||x||^2) on the chord template -void rand_exp_range(Point &lower, Point &upper, const NT &a_i, Point &p, Parameters const& var) { +void rand_exp_range(Point const &lower, Point const & upper, const NT &a_i, Point &p, Parameters const& var) { typedef typename Parameters::RNGType RNGType; NT r, r_val, fn; const NT tol = 0.00000001; @@ -314,7 +314,7 @@ void gaussian_ball_walk(Point & p, unsigned int n = P.dimension(); NT f_x, f_y, rnd; Point y = get_point_in_Dsphere(n, ball_rad); - y = y + p; + y += p; f_x = eval_exp(p, a_i); //unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); //RNGType rng(seed); diff --git a/include/samplers/samplers.h b/include/samplers/samplers.h index 97e40f186..a9f28fec0 100644 --- a/include/samplers/samplers.h +++ b/include/samplers/samplers.h @@ -44,7 +44,7 @@ Point get_direction(const unsigned int dim) { template Point get_point_on_Dsphere(const unsigned int dim, const NT &radius){ Point p = get_direction(dim); - p = (radius == 0) ? p : radius * p; + if (radius != 0) p *= radius; return p; }