From bbab7ced991d4bfd21971553093c16e904ead08f Mon Sep 17 00:00:00 2001 From: Panagiotis Repouskos Date: Mon, 23 Mar 2020 15:03:17 +0200 Subject: [PATCH] Use eigen types as internal structures (#29) * 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. * fixed bug * 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). * Eigen clean existing code, make changes in leftover files with previous implementation * Optimizations - avoid creating copies - code cleanup - more coherent coding, using Eigen and eradicating std::vector * Optimizations - Vectorize when possible the operations in line_intersect functions - Change the matrix in HPolytope to be RowMajor * 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 * Bug Error in creating Point * bug * use cooling balls in test * edit initialization * edit /test/CMakeLists.txt * 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 * leftovers from merge - make changes in code from last merge - delete tests I had previously added * cleanup - requested changes * update copyrights * bug * bug * requested changes * use += *= operators with points * use += *= operators with points --- R-proj/src/direct_sampling.cpp | 2 +- R-proj/src/sample_points.cpp | 19 +- include/annealing/gaussian_annealing.h | 7 +- include/annealing/hpoly_annealing.h | 2 +- include/annealing/ratio_estimation.h | 14 +- include/cartesian_geom/point.h | 195 ++++++++++-------- include/convex_bodies/ball.h | 46 ++--- include/convex_bodies/ballintersectconvex.h | 16 +- include/convex_bodies/ellipsoids.h | 9 +- include/convex_bodies/hpolytope.h | 177 ++++++++-------- include/convex_bodies/vpolyintersectvpoly.h | 21 +- include/convex_bodies/vpolytope.h | 66 +++--- include/convex_bodies/zonoIntersecthpoly.h | 22 +- include/convex_bodies/zpolytope.h | 32 ++- include/generators/h_polytopes_gen.h | 7 +- .../generators/known_polytope_generators.h | 1 + include/lp_oracles/solve_lp.h | 19 +- include/lp_oracles/vpolyoracles.h | 1 + include/lp_oracles/zpolyoracles.h | 1 + include/rotating.h | 1 + include/rounding.h | 21 +- include/samplers/gaussian_samplers.h | 22 +- include/samplers/samplers.h | 92 +++++---- include/volume/cooling_balls.h | 7 +- include/volume/cooling_hpoly.h | 1 + include/volume/exact_vols.h | 2 + include/volume/volume.h | 11 +- test/CMakeLists.txt | 11 +- 28 files changed, 445 insertions(+), 380 deletions(-) diff --git a/R-proj/src/direct_sampling.cpp b/R-proj/src/direct_sampling.cpp index 195cb5467..2c947b14c 100644 --- a/R-proj/src/direct_sampling.cpp +++ b/R-proj/src/direct_sampling.cpp @@ -117,7 +117,7 @@ Rcpp::NumericMatrix direct_sampling(Rcpp::Nullable body = R_NilValue unsigned int jj = 0; for (typename std::list::iterator rpit = randPoints.begin(); rpit!=randPoints.end(); rpit++, jj++) { - RetMat.col(jj) = Eigen::Map(&(*rpit).get_coeffs()[0], (*rpit).dimension()); + RetMat.col(jj) = rpit->getCoefficients(); } return Rcpp::wrap(RetMat); } diff --git a/R-proj/src/sample_points.cpp b/R-proj/src/sample_points.cpp index 09e4d3bc3..4194bd709 100644 --- a/R-proj/src/sample_points.cpp +++ b/R-proj/src/sample_points.cpp @@ -197,6 +197,7 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P = R_NilValue } } + if (Rcpp::as(random_walk).containsElementNamed("nburns")) { nburns = Rcpp::as(Rcpp::as(random_walk)["nburns"]); if (nburns < 0) { @@ -233,7 +234,8 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P = R_NilValue HP.normalize(); if (gaussian) { StartingPoint = StartingPoint - mode; - HP.shift(Eigen::Map(&mode.get_coeffs()[0], mode.dimension())); + HP.shift(mode.getCoefficients()); + } break; } @@ -252,7 +254,7 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P = R_NilValue if (billiard && diam < 0.0) VP.comp_diam(diam); if (gaussian) { StartingPoint = StartingPoint - mode; - VP.shift(Eigen::Map(&mode.get_coeffs()[0], mode.dimension())); + VP.shift(mode.getCoefficients()); } break; } @@ -271,7 +273,7 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P = R_NilValue if (billiard && diam < 0.0) ZP.comp_diam(diam); if (gaussian) { StartingPoint = StartingPoint - mode; - ZP.shift(Eigen::Map(&mode.get_coeffs()[0], mode.dimension())); + ZP.shift(mode.getCoefficients()); } break; } @@ -296,7 +298,7 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P = R_NilValue } if (gaussian) { StartingPoint = StartingPoint - mode; - VPcVP.shift(Eigen::Map(&mode.get_coeffs()[0], mode.dimension())); + VPcVP.shift(mode.getCoefficients()); } break; } @@ -342,14 +344,17 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P = R_NilValue MT RetMat(dim, numpoints); unsigned int jj = 0; + 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(&mode.get_coeffs()[0], mode.dimension()); + + RetMat.col(jj) = rpit->getCoefficients() + mode.getCoefficients(); + } else { - RetMat.col(jj) = Eigen::Map(&(*rpit).get_coeffs()[0], (*rpit).dimension()); + RetMat.col(jj) = (*rpit).getCoefficients(); } } + return Rcpp::wrap(RetMat); } 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/hpoly_annealing.h b/include/annealing/hpoly_annealing.h index 0a9f929bb..d262d84ee 100644 --- a/include/annealing/hpoly_annealing.h +++ b/include/annealing/hpoly_annealing.h @@ -37,7 +37,7 @@ void comp_diam_hpoly_zono_inter(ZonoHP &ZHP, const MT &G, const MT &AG, const VT typename std::list::iterator rpit=randPoints.begin(); NT max_norm = 0.0, iter_norm; for ( ; rpit!=randPoints.end(); rpit++) { - iter_norm = (G*Eigen::Map(&(*rpit).get_coeffs()[0], (*rpit).dimension())).norm(); + iter_norm = (G*rpit->getCoefficients()).norm(); if (iter_norm > max_norm) max_norm = iter_norm; } diams_inter.push_back(2.0 * max_norm); diff --git a/include/annealing/ratio_estimation.h b/include/annealing/ratio_estimation.h index efe5577c0..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 @@ -28,7 +29,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 +96,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 = "< e * std::abs(*mit) || - std::abs(*mit - *pit) > e * std::abs(*pit)) + for (i=0 ; icoeffs(i) - coeffs(i)) > e *std::abs(this->coeffs(i)) || + std::abs(this->coeffs(i) - coeffs(i)) > e *std::abs(coeffs(i))) return false; } return true; } + FT dot(const point& p) const { + return coeffs.dot(p.getCoefficients()); + } - FT dot(point& p){ - FT res=FT(0.0); - - typename Coeff::iterator pit=p.iter_begin(); - typename Coeff::iterator mit=coeffs.begin(); - for( ; pitcoeffs.dot(coeffs); } - - - FT squared_length() { + + + FT squared_length() const { FT lsq = FT(0.0); - typename Coeff::iterator mit=coeffs.begin(); - for ( ; mit != coeffs.end(); mit++){ - lsq += (*mit)*(*mit); + for (int i=0; i -point operator* (typename K::FT const& k, point& p) { +point operator* (const typename K::FT& k, point& p) { return p * k; } #endif + + diff --git a/include/convex_bodies/ball.h b/include/convex_bodies/ball.h index 1a470e6bb..91a7cfb82 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 @@ -15,6 +17,7 @@ struct Ball{ typedef Point BallPoint; typedef typename Point::FT NT; typedef typename std::vector::iterator viterator; + typedef Eigen::Matrix VT; Ball() {} @@ -44,26 +47,22 @@ struct Ball{ std::pair line_intersect(Point &r, Point &v) { - viterator vit=v.iter_begin(), cit=c.iter_begin(), rcit=r.iter_begin(); NT vrc(0), v2(0), rc2(0); - for( ; cit < c.iter_end() ; ++cit, ++rcit, ++vit){ - vrc += *vit * (*rcit); - v2 += *vit * (*vit); - rc2 += *rcit * (*rcit); - } + + vrc = v.dot(r); + v2 = v.dot(v); + rc2 = r.dot(r); NT disc_sqrt = std::sqrt(std::pow(vrc,2) - v2 * (rc2 - R)); return std::pair ((NT(-1)*vrc + disc_sqrt)/v2, (NT(-1)*vrc - disc_sqrt)/v2); } - std::pair line_intersect(Point &r, Point &v, const std::vector &Ar, - const std::vector &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 std::vector &Ar, - const std::vector &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); } @@ -71,24 +70,22 @@ struct Ball{ return std::pair(line_intersect(r, v).first, 0); } - std::pair line_positive_intersect(Point &r, Point &v, const std::vector &Ar, - const std::vector &Av){ + std::pair line_positive_intersect(Point &r, Point &v, const VT &Ar, + const VT &Av){ return line_positive_intersect(r, v); } - std::pair line_positive_intersect(Point &r, Point &v, const std::vector &Ar, - const std::vector &Av, NT &lambda_prev){ + std::pair line_positive_intersect(Point &r, Point &v, const VT &Ar, + const VT &Av, NT &lambda_prev){ return line_positive_intersect(r, v); } std::pair line_intersect_coord(Point &r, const unsigned int &rand_coord) { - viterator rcit=r.iter_begin(); - NT vrc = *(rcit + rand_coord); + NT vrc = r[rand_coord]; NT rc2(R); - for( ; rcit < r.iter_end() ; ++rcit){ - rc2 -= *rcit * (*rcit); - } + rc2 -= r.dot(r); + NT disc_sqrt = std::sqrt(std::pow(vrc,2) + rc2); return std::pair (NT(-1)*vrc + disc_sqrt, NT(-1)*vrc - disc_sqrt); @@ -96,7 +93,7 @@ struct Ball{ } std::pair line_intersect_coord(Point &r, const unsigned int &rand_coord, - const std::vector &lamdas) { + const VT &lamdas) { return line_intersect_coord(r, rand_coord); } @@ -104,7 +101,7 @@ struct Ball{ 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); } @@ -115,10 +112,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/ballintersectconvex.h b/include/convex_bodies/ballintersectconvex.h index 6508c9562..132a9a0c0 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 @@ -19,6 +20,7 @@ class BallIntersectPolytope { public: typedef typename CBall::NT NT; typedef typename CBall::BallPoint Point; + typedef Eigen::Matrix VT; BallIntersectPolytope() {} @@ -57,15 +59,15 @@ class BallIntersectPolytope { std::max(polypair.second, ballpair.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 = P.line_intersect(r, v, Ar, Av); std::pair ballpair = B.line_intersect(r, v); return std::pair(std::min(polypair.first, ballpair.first), std::max(polypair.second, ballpair.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 = P.line_intersect(r, v, Ar, Av, lambda_prev); std::pair ballpair = B.line_intersect(r, v); @@ -73,7 +75,7 @@ class BallIntersectPolytope { std::max(polypair.second, ballpair.second)); } - 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 = P.line_positive_intersect(r, v, Ar, Av); std::pair ball_lambda = B.line_positive_intersect(r, v); @@ -85,7 +87,7 @@ class BallIntersectPolytope { } - 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, NT &lambda_prev) { std::pair polypair = P.line_positive_intersect(r, v, Ar, Av, lambda_prev); @@ -100,7 +102,7 @@ class BallIntersectPolytope { //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 = P.line_intersect_coord(r, rand_coord, lamdas); std::pair ballpair = B.line_intersect_coord(r, rand_coord); @@ -113,7 +115,7 @@ class BallIntersectPolytope { const Point &r_prev, const unsigned int &rand_coord, const unsigned int &rand_coord_prev, - std::vector &lamdas) { + VT &lamdas) { std::pair polypair = P.line_intersect_coord(r, r_prev, rand_coord, rand_coord_prev, lamdas); std::pair ballpair = B.line_intersect_coord(r, rand_coord); diff --git a/include/convex_bodies/ellipsoids.h b/include/convex_bodies/ellipsoids.h index 21926003b..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 @@ -38,13 +39,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/hpolytope.h b/include/convex_bodies/hpolytope.h index e883a6e10..0c7e670cf 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 @@ -252,12 +253,14 @@ class HPolytope{ int is_in(const Point &p) const { 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]; + const NT* b_data = b.data(); + for (int i = 0; i < m; i++) { //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; } @@ -277,30 +280,31 @@ 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.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(); 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); - } - 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; } + + sum_nom_data++; + sum_denom_data++; } return std::pair(min_plus, max_minus); } @@ -308,73 +312,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); - NT sum_nom, sum_denom, mult; + 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(); - - for (int i = 0; i < m; i++, ++Ariter, ++Aviter) { - sum_nom = NT(0); - 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); - } - (*Ariter) = -sum_nom; - (*Aviter) = sum_denom; - sum_nom += b(i); - if (sum_denom == NT(0)) { + + Ar.noalias() = A * r.getCoefficients(); + sum_nom = b - Ar; + Av.noalias() = A * v.getCoefficients();; + + + NT* Av_data = Av.data(); + NT* sum_nom_data = sum_nom.data(); + + 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); - NT sum_nom, sum_denom, 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(); - - for (int i = 0; i < m; i++, ++Ariter, ++Aviter) { - (*Ariter) += lambda_prev * (*Aviter); - sum_nom = b(i) - (*Ariter); - sum_denom = NT(0); - j = 0; - vit = v.iter_begin(); - for ( ; vit != v.iter_end(); vit++, j++) sum_denom += A(i, j) * (*vit); - - (*Aviter) = sum_denom; - if (sum_denom == NT(0)) { + + Ar.noalias() += lambda_prev*Av; + sum_nom = b - Ar; + Av.noalias() = A * v.getCoefficients(); + + NT* sum_nom_data = sum_nom.data(); + NT* Av_data = Av.data(); + + 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); @@ -383,48 +387,48 @@ 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); - NT sum_nom, sum_denom; + 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(); + + NT* lamda_data = lamdas.data(); + NT* sum_denom_data = sum_denom.data(); for (int i = 0; i < m; i++) { - sum_nom = b(i); - sum_denom = A(i, rand_coord); - rit = r.iter_begin(); - j = 0; - for (; rit != r.iter_end(); rit++, j++) { - sum_nom -= A(i, j) * (*rit); - } - lamdas[i] = sum_nom; - 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); } @@ -435,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.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++) { - 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); } @@ -471,7 +474,7 @@ class HPolytope{ // shift polytope by a point c void shift(const VT &c){ - b = b - A*c; + b -= A*c; } @@ -510,11 +513,11 @@ 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())); - s = ((-2.0 * v.dot(s)) * s); - v = s + v; - +// 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 += -2 * v.dot(A.row(facet)) * A.row(facet); } void free_them_all() {} diff --git a/include/convex_bodies/vpolyintersectvpoly.h b/include/convex_bodies/vpolyintersectvpoly.h index e49c0d45d..c1bb9a3eb 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 @@ -210,16 +211,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 +235,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 +251,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 +265,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 cce85a7fd..14d0135f8 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 @@ -165,15 +166,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; } @@ -222,17 +219,15 @@ class VPolytope{ std::pair result; for (j=1; jgetCoefficients() - p0.getCoefficients(); } Bg = B; B = B.inverse(); for (i=0; i Ap(_d,randPoints.size()); typename std::list::iterator rpit=randPoints.begin(); - typename std::vector::iterator qit; + unsigned int i, j = 0; for ( ; rpit!=randPoints.end(); rpit++, j++) { - qit = (*rpit).iter_begin(); i=0; - for ( ; qit!=(*rpit).iter_end(); qit++, i++){ - Ap(i,j)=double(*qit); + const NT* point_data = rpit->getCoefficients().data(); + + for ( i=0; i < rpit->dimension(); i++){ + Ap(i,j)=double(*point_data); + point_data++; } } boost::numeric::ublas::matrix Q(_d, _d); @@ -355,16 +352,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); } @@ -374,14 +371,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); } @@ -390,7 +387,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; @@ -405,7 +402,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); } @@ -441,14 +438,10 @@ class VPolytope{ return false; } unsigned int j; - std::vector temp(_d,NT(0)); + typename std::vector::iterator pointIt; for (int i=0; i 1.0) a = -a; - a = a/a.norm(); - - Point s(_d, std::vector(&a[0], a.data()+a.cols()*a.rows())); + a /= a.norm(); - 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/zonoIntersecthpoly.h b/include/convex_bodies/zonoIntersecthpoly.h index c39182efb..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 @@ -63,16 +65,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 +83,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 +92,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 +104,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 +120,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 7997eafdb..e79f5a122 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 @@ -291,27 +292,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); } @@ -320,7 +321,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; @@ -337,7 +338,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); } @@ -386,17 +387,14 @@ class Zonotope { } VT a = Fmat.fullPivLu().kernel(); - NT sum = 0.0; - for (int k = 0; k < _d; ++k) sum += a(k)*p[k]; - if(sum<0.0) a = -1.0*a; + if(p.getCoefficients().dot(a) < 0.0) a = -1.0*a; a = a/a.norm(); - Point s(_d, std::vector(&a[0], a.data()+a.cols()*a.rows())); - - 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/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/generators/known_polytope_generators.h b/include/generators/known_polytope_generators.h index df52215ba..f68060238 100644 --- a/include/generators/known_polytope_generators.h +++ b/include/generators/known_polytope_generators.h @@ -306,6 +306,7 @@ Polytope gen_skinny_cube(const unsigned int &dim, bool Vpoly = false) { return P; } + /* * ToDo: brkhoff polytope generator template diff --git a/include/lp_oracles/solve_lp.h b/include/lp_oracles/solve_lp.h index 094e777e2..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 @@ -242,11 +243,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 - typename std::vector::iterator pit = direction.iter_begin(); - for(j=0; j #include #include diff --git a/include/lp_oracles/zpolyoracles.h b/include/lp_oracles/zpolyoracles.h index 3b8e853ef..06ad2c686 100644 --- a/include/lp_oracles/zpolyoracles.h +++ b/include/lp_oracles/zpolyoracles.h @@ -1,6 +1,7 @@ #ifndef ZPOLYORACLES_H #define ZPOLYORACLES_H + #include #include #include diff --git a/include/rotating.h b/include/rotating.h index a29c705cb..366b8c5f1 100644 --- a/include/rotating.h +++ b/include/rotating.h @@ -3,6 +3,7 @@ // Copyright (c) 20012-2019 Vissarion Fisikopoulos // Copyright (c) 2019 Apostolos Chalkis + // Licensed under GNU LGPL.3, see LICENCE file #ifndef ROTATING_H diff --git a/include/rounding.h b/include/rounding.h index 90e5b8f5d..63cd369f6 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 @@ -26,6 +27,7 @@ std::pair rounding_min_ellipsoid(Polytope &P , const std::pair rounding_min_ellipsoid(Polytope &P , const std::pair 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); + 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 @@ -106,18 +107,16 @@ 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(); - 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); + for (int j=0 ; rpit!=randPoints.end(); rpit++, j++) { + for (int i=0 ; idimension(); i++){ + Ap(i,j)=double((*rpit)[i]); } } boost::numeric::ublas::matrix Q(n,n); diff --git a/include/samplers/gaussian_samplers.h b/include/samplers/gaussian_samplers.h index 25da7455a..607e240fe 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 @@ -59,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; @@ -142,14 +143,14 @@ NT rand_exp_range_coord(const NT &l, const NT &u, const NT &a_i, Parameters cons // compute the first coordinate point -template +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 +176,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 +224,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 +286,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); @@ -310,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 09ce2cb8d..a9f28fec0 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 @@ -11,26 +12,30 @@ #define RANDOM_SAMPLERS_H + + // Pick a random direction as a normilized vector 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 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; } @@ -55,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; } @@ -67,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; } @@ -111,12 +116,18 @@ void boundary_rand_point_generator(Polytope &P, { typedef typename Parameters::RNGType RNGType; typedef typename Point::FT NT; + typedef typename Point::Coeff 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)); + 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); @@ -132,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); @@ -158,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); } } @@ -184,7 +195,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 typename Point::Coeff 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); @@ -201,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); @@ -221,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); @@ -250,7 +266,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 typename Point::Coeff 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); @@ -268,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); @@ -287,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); @@ -300,14 +322,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; @@ -332,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); } @@ -355,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); @@ -364,14 +386,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; @@ -391,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); @@ -427,22 +449,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; @@ -463,7 +485,7 @@ void billiard_walk(ConvexBody &P, Point &p, NT diameter, std::vector &Ar, st 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); } @@ -478,7 +500,7 @@ void billiard_walk(ConvexBody &P, Point &p, NT diameter, std::vector &Ar, st } 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++; diff --git a/include/volume/cooling_balls.h b/include/volume/cooling_balls.h index 24436ceb8..399648070 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 @@ -65,8 +67,9 @@ NT vol_cooling_balls(Polytope &P, UParameters &var, AParameters &var_ban, std::p // Save the radius of the Chebychev ball 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()); - P.shift(c_e); + + P.shift(c.getCoefficients()); + if ( !get_sequence_of_polyballs(P, BallSet, ratios, N * nu, nu, lb, ub, radius, alpha, var, rmax) ){ return -1.0; diff --git a/include/volume/cooling_hpoly.h b/include/volume/cooling_hpoly.h index 57777bd5a..7b0025900 100644 --- a/include/volume/cooling_hpoly.h +++ b/include/volume/cooling_hpoly.h @@ -3,6 +3,7 @@ // Copyright (c) 20012-2018 Vissarion Fisikopoulos // Copyright (c) 2018 Apostolos Chalkis + #ifndef COOLING_HPOLY_H #define COOLING_HPOLY_H 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 f5aeac3a3..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 @@ -78,8 +79,7 @@ NT volume(Polytope &P, } // 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()); - P.shift(c_e); + P.shift(c.getCoefficients()); c=Point(n); rnum=rnum/n_threads; @@ -288,8 +288,7 @@ NT volume_gaussian_annealing(Polytope &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()); - P.shift(c_e); + P.shift(c.getCoefficients()); // Initialization for the schedule annealing std::vector a_vals; @@ -321,7 +320,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 8478bf30e..d93d0d185 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -33,7 +33,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) @@ -50,7 +50,13 @@ else () include_directories (BEFORE ../include/lp_oracles) include_directories (BEFORE ../include/misc) - + #for Eigen + if (${CMAKE_VERSION} VERSION_LESS "3.12.0") + add_compile_options(-D "EIGEN_NO_DEBUG") + else () + add_compile_definitions("EIGEN_NO_DEBUG") + endif () + 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 +73,7 @@ else () add_library(test_main OBJECT test_main.cpp) add_executable (volume_test volume_test.cpp $) + add_executable (cheb_test chebychev_test.cpp $) #add_executable (rounding_test rounding_test.cpp $) add_executable (volumeCG_test volumeCG_test.cpp $)