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

@@ -117,7 +117,7 @@ Rcpp::NumericMatrix direct_sampling(Rcpp::Nullable<Rcpp::List> body = R_NilValue
     unsigned int jj = 0;
 
     for (typename std::list<Point>::iterator rpit = randPoints.begin(); rpit!=randPoints.end(); rpit++, jj++) {
-        RetMat.col(jj) = Eigen::Map<VT>(&(*rpit).get_coeffs()[0], (*rpit).dimension());
+        RetMat.col(jj) = rpit->getCoefficients();
     }
     return Rcpp::wrap(RetMat);
 }

R-proj/src/sample_points.cpp

@@ -197,6 +197,7 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable<Rcpp::Reference> P = R_NilValue
         }
     }
 
+
     if (Rcpp::as<Rcpp::List>(random_walk).containsElementNamed("nburns")) {
         nburns = Rcpp::as<int>(Rcpp::as<Rcpp::List>(random_walk)["nburns"]);
         if (nburns < 0) {
@@ -233,7 +234,8 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable<Rcpp::Reference> P = R_NilValue
             HP.normalize();
             if (gaussian) {
                 StartingPoint = StartingPoint - mode;
-                HP.shift(Eigen::Map<VT>(&mode.get_coeffs()[0], mode.dimension()));
+                HP.shift(mode.getCoefficients());
+
             }
             break;
         }
@@ -252,7 +254,7 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable<Rcpp::Reference> P = R_NilValue
             if (billiard && diam < 0.0) VP.comp_diam(diam);
             if (gaussian) {
                 StartingPoint = StartingPoint - mode;
-                VP.shift(Eigen::Map<VT>(&mode.get_coeffs()[0], mode.dimension()));
+                VP.shift(mode.getCoefficients());
             }
             break;
         }
@@ -271,7 +273,7 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable<Rcpp::Reference> P = R_NilValue
             if (billiard && diam < 0.0) ZP.comp_diam(diam);
             if (gaussian) {
                 StartingPoint = StartingPoint - mode;
-                ZP.shift(Eigen::Map<VT>(&mode.get_coeffs()[0], mode.dimension()));
+                ZP.shift(mode.getCoefficients());
             }
             break;
         }
@@ -296,7 +298,7 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable<Rcpp::Reference> P = R_NilValue
             }
             if (gaussian) {
                 StartingPoint = StartingPoint - mode;
-                VPcVP.shift(Eigen::Map<VT>(&mode.get_coeffs()[0], mode.dimension()));
+                VPcVP.shift(mode.getCoefficients());
             }
             break;
         }
@@ -342,14 +344,17 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable<Rcpp::Reference> P = R_NilValue
     MT RetMat(dim, numpoints);
     unsigned int jj = 0;
 
+
     for (typename std::list<Point>::iterator rpit = randPoints.begin(); rpit!=randPoints.end(); rpit++, jj++) {
         if (gaussian) {
-            RetMat.col(jj) = Eigen::Map<VT>(&(*rpit).get_coeffs()[0], (*rpit).dimension()) +
-                             Eigen::Map<VT>(&mode.get_coeffs()[0], mode.dimension());
+
+            RetMat.col(jj) = rpit->getCoefficients() + mode.getCoefficients();
+
         } else {
-            RetMat.col(jj) = Eigen::Map<VT>(&(*rpit).get_coeffs()[0], (*rpit).dimension());
+            RetMat.col(jj) = (*rpit).getCoefficients();
         }
     }
+
     return Rcpp::wrap(RetMat);
 
 }

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<NT> lamdas(P.num_of_hyperplanes(), NT(0));
+    typedef Eigen::Matrix<NT,Eigen::Dynamic,1> 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){

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<Point>::iterator rpit=randPoints.begin();
     NT max_norm = 0.0, iter_norm;
     for ( ; rpit!=randPoints.end(); rpit++) {
-        iter_norm = (G*Eigen::Map<VT>(&(*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);

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<NT>::lowest(), max_val = std::numeric_limits<NT>::max(), val, lambda;
     size_t totCount = Ntot, countIn = Ntot * ratio;
-    std::vector<NT> last_W(W), lamdas(Pb1.num_of_hyperplanes()), Av(Pb1.num_of_hyperplanes());
+    std::vector<NT> last_W(W);
+    typedef Eigen::Matrix<NT,Eigen::Dynamic,1> VT;
+    VT lamdas, Av;
+    lamdas.setZero(Pb1.num_of_hyperplanes());
+    Av.setZero(Pb1.num_of_hyperplanes());
+
     std::list<Point> randPoints;
     typename std::vector<NT>::iterator minmaxIt;
     typename std::list<Point>::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<NT> last_W(W), lamdas(Pb1.num_of_hyperplanes()), Av(Pb1.num_of_hyperplanes());
+    std::vector<NT> last_W(W);
+    typedef Eigen::Matrix<NT,Eigen::Dynamic,1> 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 = "<<countIn<<", totCount = "<<totCount<<std::endl;