Eigen implementation

R-proj/src/sample_points.cpp

@@ -339,7 +339,7 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable<Rcpp::Reference> P = 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);
 
 }

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

@@ -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<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 +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<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;

include/cartesian_geom/point.h

@@ -10,161 +10,186 @@
 #define POINT_H
 
 #include <iostream>
+#include "../external/Eigen/Eigen"
 
-template <typename K>
+template <class K>
 class point
 {
 private:
     unsigned int d;
-    typedef std::vector<typename K::FT> Coeff;
+
+    typedef Eigen::Matrix<typename K::FT, Eigen::Dynamic,1> Coeff;
+
     Coeff coeffs;
     typedef typename std::vector<typename K::FT>::iterator iter;
 public:
     typedef typename K::FT 	FT;
 
     point() {}
-    
-    point(const unsigned int &dim) {
+
+    point(const unsigned int dim) {
         d = dim;
-        coeffs = Coeff(d,0);
+        coeffs.setZero(d);
     }
-    
-    point(const unsigned int &dim, iter begin, iter endit) {
+
+    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;
     }
 
-    point(const unsigned int dim, Coeff cofs) {
+    point(const Coeff& coeffs) {
+            d = coeffs.rows();
+            this->coeffs = coeffs;
+    }
+
+    point(const unsigned int dim, std::vector<typename K::FT> cofs) {
         d = dim;
-        coeffs.assign(cofs.begin(), cofs.end());
+        iter it = cofs.begin();
+        int i=0;
+
+        for (; it != cofs.end(); it++,i++)
+            coeffs(i) = *it;
+
+    }
+
+    void add(const Coeff& coeffs) {
+        this->coeffs += coeffs;
     }
-
+
+    const Coeff& getCoefficients() const {
+        return coeffs;
+    }
+
     int dimension() const {
         return d;
     }
-    
-    void set_dimension(const unsigned int &dim) {
+
+    void set_dimension(const unsigned int dim) {
         d = dim;
     }
-    
-    void set_coord(const unsigned int &i, const FT &coord) {
-        coeffs[i] = coord;
+
+    void set_coord(const unsigned int i, FT coord) {
+        coeffs(i) = coord;
     }
 
-    Coeff get_coeffs() {
-        return coeffs;
+    FT operator[] (const unsigned int i) const {
+        return coeffs(i);
     }
-
-    FT operator[] (const unsigned int &i) const {
-        return coeffs[i];
+
+    point operator+ (const point& p) const {
+        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 operator+= (const point& p) {
+        coeffs += p.getCoefficients();
+    }
 
-        for (; pit < p.iter_end(); ++pit, ++mit, ++tmit) {
-            (*tmit) = (*pit) + (*mit);
-        }
-        return temp;
+    void operator+= (const Coeff& coeffs) {
+        this->coeffs = coeffs + this->coeffs;
     }
-
-    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 operator= (const Coeff& coeffs) {
+        this->coeffs = coeffs;
+    }
 
-        for (; pit < p.iter_end(); ++pit, ++mit, ++tmit) {
-            (*tmit) = (*mit) - (*pit);
-        }
+    point operator- (const point& p) const {
+        point temp;
+        temp.d = d;
+        temp.coeffs = coeffs - p.getCoefficients();
         return temp;
     }
 
-    point operator* (const FT& k) {
-        point temp(d, iter_begin(), iter_end());
+    point operator* (const FT k) const {
+        point temp;
+        temp.d = d;
+        temp.coeffs = coeffs * k;
+        return temp;
+    }
 
-        typename Coeff::iterator tmit = temp.iter_begin();
+    void operator*= (const FT k) {
+        coeffs *= k;
+    }
 
-        for (; tmit < temp.iter_end(); ++tmit) {
-            (*tmit) *= k;
-        }
+    point operator/ (const FT k) const {
+        point temp;
+        temp.d = d;
+        temp.coeffs = coeffs / k;
         return temp;
     }
 
+    void operator/= (const FT k)  {
+        coeffs /= k;
+    }
 
-    bool operator== (point& p) {
 
-        typename Coeff::iterator pit = p.iter_begin();
-        typename Coeff::iterator mit = coeffs.begin();
 
-        /* degree of approximation in
-        "The art of computer programming" (vol II), p. 234, Donald. E. Knuth. */
+    bool operator== (point& p) const {
+        int i=0;
+        const Coeff & coeffs = p.getCoefficients();
         FT e = 0.00000000001;
-
-        for ( ;  pit!=p.iter_end(); ++pit, ++mit) {
-            if (std::abs(*mit - *pit) > e * std::abs(*mit) ||
-                std::abs(*mit - *pit) > e * std::abs(*pit))
+        for (i=0 ; i<d ; i++) {
+            if (std::abs(this->coeffs(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(point& p){
-        FT res=FT(0.0);
+    FT dot(const point& p) const {
+        return coeffs.dot(p.getCoefficients());
+    }
 
-        typename Coeff::iterator pit=p.iter_begin();
-        typename Coeff::iterator mit=coeffs.begin();
-        for( ; pit<p.iter_end(); ++pit, ++mit){
-            res+=(*mit)*(*pit);
-        }
-        return res;
+    FT dot(const Coeff& coeffs) const {
+        return this->coeffs.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<d ; i++){
+            lsq += coeffs(i) * coeffs(i);
         }
         return lsq;
     }
 
-    void print(){
+    void print() const{
         for(unsigned int i=0; i<d; i++){
-            #ifdef VOLESTI_DEBUG
-            std::cout<<coeffs[i]<<" ";
-            #endif
+#ifdef VOLESTI_DEBUG
+            std::cout<<coeffs(i)<<" ";
+#endif
         }
-        #ifdef VOLESTI_DEBUG
+#ifdef VOLESTI_DEBUG
         std::cout<<"\n";
-        #endif
-    }
-
-
-    iter iter_begin() {
-        return coeffs.begin();
-    }
-
-    iter iter_end() {
-        return coeffs.end();
+#endif
     }
-
-
+
+
+//    iter iter_begin() {
+//        return coeffs.begin();
+//    }
+//
+//    iter iter_end() {
+//        return coeffs.end();
+//    }
+
 };
 
-template<typename K>
-point<K> operator* (typename K::FT const& k, point<K>& p) {
+template<class K>
+point<K> operator* (const typename K::FT& k, point<K>& p) {
     return p * k;
 }
 
 #endif
+
+