Reformatting changes and some minor improvements

R-proj/src/copula.cpp

@@ -63,7 +63,7 @@ Rcpp::NumericMatrix copula (Rcpp::Nullable<Rcpp::NumericVector> r1,
     typedef boost::mt19937    RNGType;
     typedef Eigen::Matrix<NT,Eigen::Dynamic,Eigen::Dynamic> MT;
     typedef Eigen::Matrix<NT,Eigen::Dynamic,1> VT;
-    typedef Ellipsoid<Point, MT> CopEll;
+    typedef Ellipsoid<Point> CopEll;
     unsigned int num_slices = 100, numpoints = 500000;
 
     if (m.isNotNull()) {

examples/order-polytope-basics/order_polytope.cpp

@@ -20,26 +20,14 @@ std::vector<T> linspace(T a, T b, size_t N) {
 typedef typename Poset::RT RT;
 typedef typename Poset::RV RV;
 
-Poset read_poset_from_file(std::string filename) {
-    std::ifstream data_file;
-    data_file.open(filename);
-
-    unsigned int n;
-    data_file >> n;
-
-    RT curr_relation;
-    RV relations;
-    while(data_file >> curr_relation.first >> curr_relation.second)
-        relations.push_back(curr_relation);
-
-    return Poset(n, relations);
-}
-
 
 int main(int argc, char const *argv[]) {
     std::cout << "\nPoset operations: \n";
     //  ----------- basic poset operations -----------
-    Poset poset = read_poset_from_file(std::string(argv[1]));
+    std::string filename (argv[1]);
+    std::ifstream data_file;
+    data_file.open(filename);
+    Poset poset = read_poset_from_file(data_file);
     poset.print();
 
     std::cout << "Checking if a sample Point (linearly spaced coordinates) lies inside the poset: ";
@@ -91,11 +79,11 @@ int main(int argc, char const *argv[]) {
     ray.set_coord(0, 1.5);
     std::cout << "incident ray: ";
     ray.print();
-    
+
     OP.compute_reflection(ray, Point(), 2*OP.dimension());
     std::cout << "reflected ray: ";
     ray.print();
-    // ---------------------------------------------    
+    // ---------------------------------------------
 
     return 0;
 }

include/convex_bodies/ellipsoid.h

@@ -16,23 +16,17 @@
 
 #include <iostream>
 #include <Eigen/Eigen>
-#include <boost/math/special_functions/gamma.hpp>
+#include "volume/math_helpers.hpp"
 
 
-template <typename NT>
-NT log_gamma_function(NT x)
-{
-    if (x <= NT(100)) return std::log(tgamma(x));
-    return (std::log(x - NT(1)) + log_gamma_function(x - NT(1)));
-}
-
-
-template <class Point, class MT>
+template <class Point>
 class Ellipsoid{
-private:
+public:
+typedef Point PointType;
     typedef typename Point::FT NT;
     typedef typename std::vector<NT>::iterator viterator;
     typedef Eigen::Matrix<NT, Eigen::Dynamic, 1> VT;
+    typedef Eigen::Matrix<NT, Eigen::Dynamic, Eigen::Dynamic> MT;
 
     // representation is (x - c)' A (x - c) <= 1, center is assumed to be origin for now
     MT A;
@@ -66,7 +60,7 @@ class Ellipsoid{
     }
 
 
-    // Constructor for copula ellipsoid
+    // Constructor for copula ellipsoid only
     Ellipsoid(std::vector<std::vector<NT> >& Ain) {
         _dim = Ain.size();
         A.resize(_dim, _dim);
@@ -75,20 +69,6 @@ class Ellipsoid{
                 A(i,j) = Ain[i][j];
             }
         }
-
-        Eigen::SelfAdjointEigenSolver<MT> eigensolver(A);
-        if (eigensolver.info() != Eigen::Success) {
-            throw std::runtime_error("Eigen solver returned error!");
-        }
-
-        _eigen_values = eigensolver.eigenvalues();
-        _Eigen_Vectors = eigensolver.eigenvectors();
-
-        _eigen_values_inv = _eigen_values.array().inverse().matrix();
-        _eigen_values_inv_sqrt = _eigen_values_inv.array().sqrt().matrix();
-
-        _dim = A.rows();
-        c = Point(_dim);
     }
 
 

include/convex_bodies/orderpolytope.h

@@ -93,8 +93,20 @@ class OrderPolytope {
     }
 
 
+    Eigen::SparseMatrix<NT> get_mat() const
+    {
+        return _A.sparseView();
+    }
+
+
+    VT get_vec() const
+    {
+        return b;
+    }
+
+
     // print polytope in Ax <= b format
-    void print()
+    void print() const
     {
         std::cout << " " << _A.rows() << " " << _d << " double" << std::endl;
         for (unsigned int i = 0; i < _A.rows(); i++) {
@@ -647,6 +659,34 @@ class OrderPolytope {
     }
 
 
+    // get a point inside the order polytope,
+    // NOTE: the current implementation only works for non shifted order polytope
+    VT inner_point()
+    {
+        // get topologically sorted list of indices
+        std::vector<unsigned int> sorted_list = _poset.topologically_sorted_list();
+
+        // vector to hold n linearly spaced values between 0-1
+        std::vector<NT> lin_space_values(_d);
+        NT start = 0.05, end = 0.95;
+        NT h = (end - start)/static_cast<NT>(_d-1);
+        NT val = start;
+        for(auto x=lin_space_values.begin(); x!=lin_space_values.end(); ++x) {
+            *x = val;
+            val += h;
+        }
+
+        // final result vector
+        VT res(_d);
+        for(int i=0; i<_d; ++i) {
+            unsigned int curr_idx = sorted_list[i];
+            res(curr_idx) = lin_space_values[i];
+        }
+
+        return res;
+    }
+
+
     void normalize()
     {
         if (_normalized == true)

include/misc/misc.h

@@ -1,6 +1,7 @@
 // VolEsti (volume computation and sampling library)
 
 // Copyright (c) 2012-2018 Vissarion Fisikopoulos
+// Copyright (c) 2021 Vaibhav Thakkar
 
 // Licensed under GNU LGPL.3, see LICENCE file
 
@@ -9,6 +10,9 @@
 
 #include <iostream>
 #include <vector>
+#include <fstream>
+#include <sstream>
+#include "poset.h"
 
 //function to print rounding to double coordinates
 template <class T>
@@ -203,4 +207,62 @@ std::pair<Point, NT> read_inner_ball(std::istream &is) {
     return std::make_pair(center, radius);
 }
 
+/* read a poset given in the following format:
+    - First line contains a single positive integer 'n' - number of elements
+    - Next `m` lines follow containing a pair 'i j' in each line to signify A_i <= A_j
+        i.e i_th element is less than or equal to the j_th element
+*/
+Poset read_poset_from_file(std::istream &data_file) {
+    typedef typename Poset::RT RT;
+    typedef typename Poset::RV RV;
+
+    // read number of elements
+    unsigned int n;
+    data_file >> n;
+
+    // read relations line by line
+    RT curr_relation;
+    RV relations;
+    while(data_file >> curr_relation.first >> curr_relation.second)
+        relations.push_back(curr_relation);
+
+    return Poset(n, relations);
+}
+
+
+// read a poset given as an adjacency matrix
+std::pair<bool, Poset> read_poset_from_file_adj_matrix(std::istream &in) {
+    typedef typename Poset::RV RV;
+
+    RV edges;
+    unsigned int x, n = 0;
+
+    // read a single line
+    std::string line;
+    std::getline(in, line);
+    std::stringstream line_ss(line);
+    while(line_ss >> x) {
+        if(x) {
+            edges.emplace_back(0, n);
+        }
+        ++n;
+    }
+
+    // read rest of the lines
+    for(unsigned int a = 1; a < n; ++a) {
+        for(unsigned int b = 0; b < n; ++b) {
+            if(!(in >> x)) {
+                std::cerr << "Invalid adjacency matrix";
+                return std::pair<bool, Poset>(false, Poset());
+            }
+
+            if(x) {
+                edges.emplace_back(a, b);
+            }
+        }
+    }
+
+    return std::pair<bool, Poset>(true, Poset(n, edges));
+}
+
 #endif //MISC_H

include/misc/poset.h

@@ -13,6 +13,7 @@
 
 #include <iostream>
 #include <vector>
+#include <queue>
 
 
 class Poset {
@@ -25,6 +26,8 @@ class Poset {
     RV order_relations;     // pairs of form a <= b
 
 public:
+    Poset() {}
+
     Poset(unsigned int _n, RV& _order_relations) :
         n(_n), order_relations(verify(_order_relations, n))
     {}
@@ -89,6 +92,40 @@ class Poset {
 
         return true;
     }
+
+
+    std::vector<unsigned int> topologically_sorted_list() const
+    {
+        std::vector<std::vector<unsigned int> > adjList(n);
+        std::vector<unsigned int> indegree(n, 0);
+
+        for(auto x: order_relations) {
+            adjList[x.first].push_back(x.second);
+            indegree[x.second] += 1;
+        }
+
+        std::queue<unsigned int> q;
+        for(unsigned int i=0; i<n; ++i) {
+            if(indegree[i] == 0)
+                q.push(i);
+        }
+
+        std::vector<unsigned int> res;
+        while(!q.empty()) {
+            unsigned int curr = q.front();
+            res.push_back(curr);
+            q.pop();
+
+            for(unsigned int i=0; i<adjList[curr].size(); ++i) {
+                unsigned int adj_idx = adjList[curr][i];
+                indegree[adj_idx] -= 1;
+                if(indegree[adj_idx] == 0)
+                    q.push(adj_idx);
+            }
+        }
+
+        return res;
+    }
 };
 
 #endif

include/random_walks/uniform_billiard_walk.hpp

@@ -22,6 +22,7 @@
     #include "convex_bodies/zonoIntersecthpoly.h"
 #endif
 #include "sampling/sphere.hpp"
+#include "random_walks/boundary_cdhr_walk.hpp"
 #include "generators/boost_random_number_generator.hpp"
 #include "sampling/random_point_generators.hpp"
 #include "volume/sampling_policies.hpp"

include/volume/math_helpers.hpp

@@ -0,0 +1,46 @@
+// VolEsti (volume computation and sampling library)
+
+// Copyright (c) 2012-2020 Vissarion Fisikopoulos
+// Copyright (c) 2018-2020 Apostolos Chalkis
+// Copyright (c) 2021 Vaibhav Thakkar
+
+// Contributed and/or modified by Vaibhav Thakkar, as part of Google Summer of Code 2021 program.
+
+// Licensed under GNU LGPL.3, see LICENCE file
+
+#ifndef MATH_HELPERS_HPP
+#define MATH_HELPERS_HPP
+
+#include <vector>
+#include <boost/math/special_functions/gamma.hpp>
+
+
+//An implementation of Welford's algorithm for mean and variance.
+template <typename NT>
+std::pair<NT, NT> get_mean_variance(std::vector<NT>& vec)
+{
+    NT mean = 0;
+    NT M2 = 0;
+    NT variance = 0;
+    NT delta;
+
+    unsigned int i=0;
+    for (auto vecit = vec.begin(); vecit!=vec.end(); vecit++, i++)
+    {
+        delta = *vecit - mean;
+        mean += delta / (i + 1);
+        M2 += delta * (*vecit - mean);
+        variance = M2 / (i + 1);
+    }
+    return std::pair<NT, NT> (mean, variance);
+}
+
+
+template <typename NT>
+static NT log_gamma_function(NT x)
+{
+    if (x <= NT(100)) return std::log(tgamma(x));
+    return (std::log(x - NT(1)) + log_gamma_function(x - NT(1)));
+}
+
+#endif  // MATH_HELPERS_HPP

include/volume/volume_cooling_balls.hpp

@@ -20,8 +20,11 @@
     #include "convex_bodies/zpolytope.h"
     #include "convex_bodies/vpolyintersectvpoly.h"
 #endif
+#include "random_walks/uniform_cdhr_walk.hpp"
+#include "convex_bodies/ball.h"
 #include "convex_bodies/ballintersectconvex.h"
 #include "sampling/random_point_generators.hpp"
+#include "volume/math_helpers.hpp"
 
 
 ////////////////////////////////////
@@ -685,15 +688,6 @@ NT estimate_ratio_interval(PolyBall1 const& Pb1,
     return NT(ratio_parameters.count_in) / NT(ratio_parameters.tot_count);
 }
 
-
-template <typename NT>
-NT log_gamma_function(NT x) 
-{
-    if (x <= NT(100)) return std::log(tgamma(x));
-    return (std::log(x - NT(1)) + log_gamma_function(x - NT(1)));
-}
-
-
 template
 <
     typename WalkTypePolicy,