Skip to content

Commit

Permalink
Merge pull request #6 from TolisChal/use_eigen_matrix
Browse files Browse the repository at this point in the history
Use eigen matrix for the H-polytope representation. Able to use float or long double.
TolisChal authored Jul 31, 2018

Verified

This commit was created on GitHub.com and signed with GitHub’s verified signature. The key has expired.
2 parents 1ecad41 + de9ee55 commit 52da8e8
Showing 13 changed files with 446 additions and 611 deletions.
16 changes: 8 additions & 8 deletions R-proj/src/vol_R.cpp
Original file line number Diff line number Diff line change
@@ -17,7 +17,7 @@ double vol_R(Rcpp::NumericMatrix A, int W ,double e, Rcpp::NumericVector Chebych

int n, nexp=1, n_threads=1,i,j;
int walk_len;//to be defined after n
double exactvol(-1.0);
NT exactvol(-1.0);
bool verbose=true,
rand_only=false,
round_only=false,
@@ -30,7 +30,7 @@ double vol_R(Rcpp::NumericMatrix A, int W ,double e, Rcpp::NumericVector Chebych
rotate=false,
experiments=true,
coordinate=coord;
Polytope<double> P;
Polytope<NT> P;

walk_len=W;
n=A.ncol()-1;
@@ -52,8 +52,8 @@ double vol_R(Rcpp::NumericMatrix A, int W ,double e, Rcpp::NumericVector Chebych
boost::random::uniform_real_distribution<> urdist1(-1,1);

vars var(rnum,n,walk_len,n_threads,0.0,0.0,0,0.0,0,rng,urdist,urdist1,verbose,rand_only,round,NN,birk,coordinate);
std::vector<std::vector<double> > Pin(m+1, std::vector<double>(n+1));
std::vector<double> bin(m);
std::vector<std::vector<NT> > Pin(m+1, std::vector<NT>(n+1));
std::vector<NT> bin(m);

for (i=0; i<m+1; i++){
//bin[i]=b[i];
@@ -69,11 +69,11 @@ double vol_R(Rcpp::NumericMatrix A, int W ,double e, Rcpp::NumericVector Chebych
}
P.init(Pin);
//Compute chebychev ball//
std::pair<Point,double> CheBall;
std::pair<Point,NT> CheBall;
if(Chebychev.size()!=P.dimension()+1){
CheBall = P.chebyshev_center();
}else{
std::vector<double> temp_p;
std::vector<NT> temp_p;
for (int j=0; j<P.dimension(); j++){
temp_p.push_back(Chebychev[j]);
}
@@ -82,10 +82,10 @@ double vol_R(Rcpp::NumericMatrix A, int W ,double e, Rcpp::NumericVector Chebych
CheBall.first = xc; CheBall.second = radius;
}

Polytope<double> P_to_test(P);
//Polytope<double> P_to_test(P);

NT vol = volume(P,var,var,CheBall);


return vol;
return ((double)vol);
}
26 changes: 14 additions & 12 deletions external/LPsolve_src/run_headers/solve_lp.h
Original file line number Diff line number Diff line change
@@ -9,14 +9,14 @@
#include "lp_lib.h"


template <typename K>
std::pair<Point,double> solveLP(std::vector<std::vector<K> > A, int d){
//template <typename K>
std::pair<Point,NT> solveLP(Eigen::MatrixXd &A, Eigen::VectorXd &b, int d){

//typedef typename T1::FT K;
//int d=P.dimension();
//std::vector<std::vector<K> > A = P.get_matrix();
lprec *lp;
int Ncol=d+1, *colno = NULL, j, ret = 0, m=A.size();
int Ncol=d+1, *colno = NULL, j, ret = 0, m=A.rows();
REAL *row = NULL;

lp = make_lp(m, Ncol);
@@ -46,22 +46,24 @@ std::pair<Point,double> solveLP(std::vector<std::vector<K> > A, int d){
set_add_rowmode(lp, TRUE); /* makes building the model faster if it is done rows by row */
}
int i=0;
K sum;
NT sum;
while(ret==0 & i<m){
/* construct all rows (120 x + 210 y <= 15000) */
sum=K(0);
sum=NT(0);
for(j=0; j<d; j++){
colno[j] = j+1; /* j_th column */
row[j] = A[i][j+1];
sum += A[i][j+1]*A[i][j+1];
//row[j] = A[i][j+1];
row[j] = A(i,j);
sum+=A(i,j)*A(i,j);
//sum += A[i][j+1]*A[i][j+1];
}
colno[d] = d+1; /* last column */
row[d] = std::sqrt(sum);
//set_bounds(lp, d, 0.0, infinite);


/* add the row to lpsolve */
if(!add_constraintex(lp, d+1, row, colno, LE, A[i][0])){
if(!add_constraintex(lp, d+1, row, colno, LE, b(i))){
ret = 3;
}
i++;
@@ -107,18 +109,18 @@ std::pair<Point,double> solveLP(std::vector<std::vector<K> > A, int d){
}

//if(ret == 0) {
std::vector<double> temp_p(d,0);
std::vector<NT> temp_p(d,0);
get_variables(lp, row);
for(j = 0; j < d; j++){
//printf("%s: %f\n", get_col_name(lp, j + 1), row[j]);
temp_p[j]=double(row[j]);
temp_p[j]=NT(row[j]);
}

Point xc( d , temp_p.begin() , temp_p.end() );
double r=double(get_objective(lp));
NT r=NT(get_objective(lp));
delete_lp(lp);

return std::pair<Point,double> (xc,r);
return std::pair<Point,NT> (xc,r);


}
111 changes: 50 additions & 61 deletions include/cartesian_geom/point.h
Original file line number Diff line number Diff line change
@@ -23,107 +23,96 @@ class point

point() {}

point(const int dim){
d=dim;
coeffs=Coeff(d,0);
point(const int dim) {
d = dim;
coeffs = Coeff(d,0);
}

point(const int dim, iter begin, iter end){
d=dim;
coeffs=Coeff(begin,end);
point(const int dim, iter begin, iter end) {
d = dim;
coeffs = Coeff(begin,end);
}

int dimension(){
int dimension() {
return d;
}

void set_dimension(const int dim){
d=dim;
void set_dimension(const int dim) {
d = dim;
}

void set_coord(const int i, FT coord){
coeffs[i]=coord;
void set_coord(const int i, FT coord) {
coeffs[i] = coord;
}

FT operator[] (const int i){
FT operator[] (const int i) {
return coeffs[i];
}

point operator+ (point& p){
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);
}
//for (int i=0; i<d; i++){
//temp.coeffs[i]=p[i]+coeffs[i];
//}

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);
}
return temp;
}

point operator- (point& p){
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)=(*mit)-(*pit);
}
//for (int i=0; i<d; i++){
//temp.coeffs[i]=coeffs[i]-p[i];
//}


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) = (*mit) - (*pit);
}
return temp;
}

point operator* (const FT& k){
point temp(d,iter_begin(),iter_end());
//point temp=*this;
typename Coeff::iterator tmit=temp.iter_begin();
for( ; tmit<temp.iter_end(); ++tmit){
(*tmit)*=k;
}
//for (int i=0; i<d; i++){
// temp.coeffs[i]=coeffs[i]*k;
//}

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;
}
return temp;
}


FT squared_length(){

FT lsq=FT(0);

for (int i=0; i<d; i++){
lsq+=coeffs[i]*coeffs[i];
}
FT squared_length() {

FT lsq = FT(0);

typename Coeff::iterator mit=coeffs.begin();
for ( ; mit != coeffs.end(); mit++){
lsq += (*mit)*(*mit);
}
return lsq;
}


iter iter_begin(){

iter iter_begin() {
return coeffs.begin();

}

iter iter_end(){

iter iter_end() {
return coeffs.end();

}


};

template<class K>
point<K> operator* (const typename K::FT& k, point<K>& p){

return p*k;

point<K> operator* (const typename K::FT& k, point<K>& p) {
return p * k;
}

#endif
89 changes: 27 additions & 62 deletions include/convex_bodies/ballintersectconvex.h
Original file line number Diff line number Diff line change
@@ -93,7 +93,7 @@ struct Ball{
};


template <class T>
template <class T, typename FT>
class BallIntersectPolytope {
public:
BallIntersectPolytope(T &P, Ball &B) : _P(P), _B(B) {};
@@ -102,7 +102,6 @@ class BallIntersectPolytope {
Ball second() { return _B; }

int is_in(Point p){
//std::cout << "calling is in"<<std::endl;
if(_B.is_in(p)==-1)
return _P.is_in(p);
return 0;
@@ -116,68 +115,44 @@ class BallIntersectPolytope {
return _P.dimension();
}

std::pair<NT,NT> line_intersect(Point r,
std::pair<FT,FT> line_intersect(Point r,
Point v) {

std::pair <NT, NT> polypair = _P.line_intersect(r, v);
std::pair <NT, NT> ballpair = _B.line_intersect(r, v);
return std::pair<NT, NT>(std::min(polypair.first, ballpair.first),
std::pair <FT, FT> polypair = _P.line_intersect(r, v);
std::pair <FT, FT> ballpair = _B.line_intersect(r, v);
return std::pair<FT, FT>(std::min(polypair.first, ballpair.first),
std::max(polypair.second, ballpair.second));
}
/*
std::pair<Point,Point> returnpair;
std::pair<Point,Point> ballpair;
bool ballinter=false;

//check the first intersection point if it is inside ball
if(_B.is_in(polypair.first)){
//std::cout<<"inside ball 1, radius:"<<_B.radius()<<std::endl;
//std::cout<<polypair.first<<std::endl;
returnpair.first = polypair.first;
}else{
//std::cout<<"outside ball 1, radius:"<<_B.radius()<<std::endl;
//std::cout<<polypair.first<<std::endl;
ballinter=true;
//compute the intersection with ball
ballpair = _B.line_intersect(r,v);
returnpair.first = ballpair.first;
//std::cout<<returnpair.first<<std::endl;
}
//check the second intersection point
if(_B.is_in(polypair.second)){
//std::cout<<"inside ball 2, radius:"<<_B.radius()<<std::endl;
//std::cout<<polypair.second<<std::endl;
returnpair.second = polypair.second;
}else{
//std::cout<<"outside ball 2, radius:"<<_B.radius()<<std::endl;
//std::cout<<polypair.second<<std::endl;
if(ballinter) //if the intersection with ball is already computed
returnpair.second = ballpair.second;
else returnpair.second = (_B.line_intersect(r,v)).second;
//std::cout<<returnpair.second<<std::endl;
}
return returnpair;
}*/
//First coordinate ray shooting intersecting convex body
std::pair<FT,FT> line_intersect_coord(Point &r,
int rand_coord,
std::vector<FT> &lamdas) {

std::pair <FT, FT> polypair = _P.line_intersect_coord(r, rand_coord, lamdas);
std::pair <FT, FT> ballpair = _B.line_intersect_coord(r, rand_coord);
return std::pair<FT, FT>(std::min(polypair.first, ballpair.first),
std::max(polypair.second, ballpair.second));
}

std::pair<NT,NT> line_intersect_coord(Point &r,
//Not the first coordinate ray shooting intersecting convex body
std::pair<FT,FT> line_intersect_coord(Point &r,
Point &r_prev,
int rand_coord,
int rand_coord_prev,
std::vector<NT> &lamdas,
bool init
){
std::vector<FT> &lamdas) {

std::pair<NT,NT> polypair = _P.line_intersect_coord(r,r_prev,rand_coord,rand_coord_prev,lamdas,init);
std::pair<NT,NT> ballpair = _B.line_intersect_coord(r,rand_coord);
return std::pair<NT,NT> (std::min(polypair.first,ballpair.first),
std::max(polypair.second,ballpair.second));
std::pair <FT, FT> polypair = _P.line_intersect_coord(r, r_prev, rand_coord, rand_coord_prev, lamdas);
std::pair <FT, FT> ballpair = _B.line_intersect_coord(r, rand_coord);
return std::pair<FT, FT>(std::min(polypair.first, ballpair.first),
std::max(polypair.second, ballpair.second));
}

std::pair<NT,NT> query_dual(Point &p, int rand_coord){
std::pair<NT,NT> polypair = _P.query_dual(p,rand_coord);
std::pair<NT,NT> ballpair = _B.line_intersect_coord(p,rand_coord);
return std::pair<NT,NT> (std::min(polypair.first,ballpair.first),
std::max(polypair.second,ballpair.second));
std::pair<FT,FT> query_dual(Point &p, int rand_coord) {
std::pair <FT, FT> polypair = _P.query_dual(p, rand_coord);
std::pair <FT, FT> ballpair = _B.line_intersect_coord(p, rand_coord);
return std::pair<FT, FT>(std::min(polypair.first, ballpair.first),
std::max(polypair.second, ballpair.second));
}

private:
@@ -224,30 +199,20 @@ class PolytopeIntersectEllipsoid {

//check the first intersection point if it is inside ball
if(E.is_in(polypair.first)){
//std::cout<<"inside ball 1, radius:"<<_B.radius()<<std::endl;
//std::cout<<polypair.first<<std::endl;
returnpair.first = polypair.first;
}else{
//std::cout<<"outside ball 1, radius:"<<_B.radius()<<std::endl;
//std::cout<<polypair.first<<std::endl;
ellinter=true;
//compute the intersection with ball
ellpair = E.line_intersect(r,v);
returnpair.first = ellpair.first;
//std::cout<<returnpair.first<<std::endl;
}
//check the second intersection point
if(E.is_in(polypair.second)){
//std::cout<<"inside ball 2, radius:"<<_B.radius()<<std::endl;
//std::cout<<polypair.second<<std::endl;
returnpair.second = polypair.second;
}else{
//std::cout<<"outside ball 2, radius:"<<_B.radius()<<std::endl;
//std::cout<<polypair.second<<std::endl;
if(ellinter) //if the intersection with ball is already computed
returnpair.second = ellpair.second;
else returnpair.second = (E.line_intersect(r,v)).second;
//std::cout<<returnpair.second<<std::endl;
}
return returnpair;
}
438 changes: 188 additions & 250 deletions include/convex_bodies/polytopes.h

Large diffs are not rendered by default.

6 changes: 3 additions & 3 deletions include/misc.h
Original file line number Diff line number Diff line change
@@ -136,7 +136,7 @@ void print_polymake_volfile2(T &P,
}*/

int read_pointset(std::istream &is,
std::vector<std::vector<double> > &Input){
std::vector<std::vector<NT> > &Input){

std::string point;

@@ -152,7 +152,7 @@ int read_pointset(std::istream &is,
//ignore empty spaces on start of line
found = point.find_first_not_of(" ",found);

std::vector<double> input;
std::vector<NT> input;
while (found2!=std::string::npos || point[found]=='-')
{
//std::cout<<"*"<<(point[found]!='-')<<"*"<<std::endl;
@@ -161,7 +161,7 @@ int read_pointset(std::istream &is,
found2 = point.find_first_not_of("0123456789-",found);

//std::cout<<point.substr(found,found2-found)<<" ";
double num = atof(point.substr(found,found2-found).c_str());
NT num = atof(point.substr(found,found2-found).c_str());
found=point.find_first_not_of(" ",found2);
//std::cout<<"found"<<point[found]<<std::endl;
if(point[found]=='/'){
51 changes: 19 additions & 32 deletions include/rounding.h
Original file line number Diff line number Diff line change
@@ -28,7 +28,7 @@ Eigen::MatrixXd getPointsMat(std::list<Point> randPoints, int dim){

template <class T1>
std::pair<Point, NT> approx_R(T1 &P, vars var){
std::pair<Point,double> Cheb_ball=solveLP(P.get_matrix(), P.dimension());
std::pair<Point,double> Cheb_ball=P.chebyshev_center();
Point c=Cheb_ball.first;
NT radius = Cheb_ball.second;

@@ -156,12 +156,12 @@ NT rounding_SVD(T1 &P , Point c, NT radius, vars &var){


// ----- ROUNDING ------ //
template <class T1>
std::pair<double,double> rounding_min_ellipsoid(T1 &P , std::pair<Point,double> CheBall, vars &var){
template <class T1, typename FT>
std::pair <FT, FT> rounding_min_ellipsoid(T1 &P , std::pair<Point,FT> CheBall, vars &var) {
int n=var.n, walk_len=var.walk_steps;
bool print=var.verbose;
Point c = CheBall.first;
NT radius = CheBall.second;
FT radius = CheBall.second;
// 2. Generate the first random point in P
// Perform random walk on random point in the Chebychev ball
Random_points_on_sphere_d<Point> gen (n, radius);
@@ -174,15 +174,15 @@ std::pair<double,double> rounding_min_ellipsoid(T1 &P , std::pair<Point,double>
int num_of_samples = 10*n;//this is the number of sample points will used to compute min_ellipoid
randPoints.clear();
rand_point_generator(P, p, num_of_samples, walk_len, randPoints, var);
NT current_dist, max_dist;
FT current_dist, max_dist;
for(std::list<Point>::iterator pit=randPoints.begin(); pit!=randPoints.end(); ++pit){
current_dist=(*pit-c).squared_length();
if(current_dist>max_dist){
max_dist=current_dist;
}
}
max_dist=std::sqrt(max_dist);
NT R=max_dist/radius;
FT R=max_dist/radius;
int mm=randPoints.size();
boost::numeric::ublas::matrix<double> Ap(n,mm);
for(int j=0; j<mm; j++){
@@ -195,39 +195,31 @@ std::pair<double,double> rounding_min_ellipsoid(T1 &P , std::pair<Point,double>
boost::numeric::ublas::matrix<double> Q(n,n);
boost::numeric::ublas::vector<double> c2(n);
size_t w=1000;
double elleps=Minim::KhachiyanAlgo(Ap,0.01,w,Q,c2);
FT elleps=Minim::KhachiyanAlgo(Ap,0.01,w,Q,c2);

Eigen::MatrixXd E(n,n);
Eigen::VectorXd e(n);

//Get ellipsoid matrix and center as Eigen objects
for(int i=0; i<n; i++){
e(i)=c2(i);
e(i)=FT(c2(i));
for (int j=0; j<n; j++){
E(i,j)=Q(i,j);
E(i,j)=FT(Q(i,j));
}
}

int m=P.num_of_hyperplanes();
Eigen::MatrixXd A(m,n);
Eigen::VectorXd b(m);
for(int i=0; i<m; ++i){
b(i) = P.get_coeff(i,0);
for(int j=1; j<n+1; ++j){
A(i,j-1) = P.get_coeff(i,j);
}
}
Eigen::MatrixXd A = P.get_eigen_mat();
Eigen::VectorXd b = P.get_eigen_vec();

//Find the smallest and the largest axes of the elliposoid
Eigen::EigenSolver<Eigen::MatrixXd> eigensolver(E);
NT rel = std::real(eigensolver.eigenvalues()[0]);
NT Rel = std::real(eigensolver.eigenvalues()[0]);
FT rel = std::real(eigensolver.eigenvalues()[0]);
FT Rel = std::real(eigensolver.eigenvalues()[0]);
for(int i=1; i<n; i++){
if(std::real(eigensolver.eigenvalues()[i])<rel) rel=std::real(eigensolver.eigenvalues()[i]);
if(std::real(eigensolver.eigenvalues()[i])>Rel) Rel=std::real(eigensolver.eigenvalues()[i]);

}
//std::cout<<rel<<" "<<Rel<<std::endl;

Eigen::LLT<Eigen::MatrixXd> lltOfA(E); // compute the Cholesky decomposition of E
Eigen::MatrixXd L = lltOfA.matrixL(); // retrieve factor L in the decomposition
@@ -239,15 +231,10 @@ std::pair<double,double> rounding_min_ellipsoid(T1 &P , std::pair<Point,double>
A = A*(L_1.transpose());

// Write changes (actually perform rounding) to the polytope!
for(int i=0; i<m; ++i){
P.put_coeff(i,0,b(i));
for(int j=1; j<n+1; ++j){
P.put_coeff(i,j,A(i,j-1));
}
}
P.set_eigen_mat(A);
P.set_eigen_vec(b);

//return L_1.determinant();
return std::pair<double,double> (L_1.determinant(),rel/Rel);
return std::pair<FT,FT> (L_1.determinant(),rel/Rel);
}


@@ -311,12 +298,12 @@ double rotating_old(T &P){

// -------- ROTATION ---------- //
template <class T>
double rotating(T &P){
NT rotating(T &P){

bool print = true;
//if(print) std::cout<<"\nRotate..."<<std::endl;
typedef Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> MT;
typedef Eigen::Matrix<double,Eigen::Dynamic,1> VT;
typedef Eigen::Matrix<NT,Eigen::Dynamic,Eigen::Dynamic> MT;
typedef Eigen::Matrix<NT,Eigen::Dynamic,1> VT;

int m=P.num_of_hyperplanes();
int d=P.dimension();
209 changes: 83 additions & 126 deletions include/samplers.h
Original file line number Diff line number Diff line change
@@ -30,9 +30,9 @@ class Random_points_on_sphere_d

Random_points_on_sphere_d() {}

Random_points_on_sphere_d(int dim, FT radius){
d=dim;
r=radius;
Random_points_on_sphere_d(int dim, FT radius) {
d = dim;
r = radius;
}

template <typename GeneratorType>
@@ -54,7 +54,7 @@ class Random_points_on_sphere_d
return point;
}


};

// WARNING: USE ONLY WITH BIRKHOFF POLYOPES
@@ -69,7 +69,7 @@ int birk_sym(T &P,K &randPoints,Point &p){

//std::cout << "The n! possible permutations with n elements:\n";
do {
std::vector<double> newpv;
std::vector<NT> newpv;
for (int i=0; i<n; i++){
//std::cout << myints[i] << " ";
}
@@ -103,105 +103,91 @@ int rand_point_generator(T &P,
int rnum,
int walk_len,
K &randPoints,
vars &var // constans for volume
)
vars &var) // constans for volume
{
//std::cout<<"EDW2!"<<std::endl;
typedef typename Point::FT FT;
int n = var.n;
bool birk = var.birk;
RNGType &rng = var.rng;
boost::random::uniform_real_distribution<> urdist(0,1);
boost::random::uniform_int_distribution<> uidist(0,n-1);
//std::uniform_real_distribution<NT> urdist = var.urdist;
//std::uniform_int_distribution<int> uidist(0,n-1);

std::vector<NT> lamdas(P.num_of_hyperplanes(),NT(0));
//int rand_coord = rand()%n;
//double kapa = double(rand())/double(RAND_MAX);
int rand_coord = uidist(rng);
double kapa = urdist(rng);
boost::random::uniform_real_distribution<> urdist(0, 1);
boost::random::uniform_int_distribution<> uidist(0, n - 1);

std::vector <NT> lamdas(P.num_of_hyperplanes(), NT(0));
int rand_coord, rand_coord_prev;
FT kapa;
Point p_prev = p;
if(var.coordinate){
//std::cout<<"[1a]P dim: "<<P.dimension()<<std::endl;
hit_and_run_coord_update(p,p_prev,P,rand_coord,rand_coord,kapa,lamdas,var,var,true);
//std::cout<<"[1b]P dim: "<<P.dimension()<<std::endl;
}else
hit_and_run(p,P,var,var);

for(int i=1; i<=rnum; ++i){

for(int j=0; j<walk_len; ++j){
int rand_coord_prev = rand_coord;
//rand_coord = rand()%n;
//kapa = double(rand())/double(RAND_MAX);
rand_coord = uidist(rng);
kapa = urdist(rng);
if(var.coordinate){
//std::cout<<"[1c]P dim: "<<P.dimension()<<std::endl;
hit_and_run_coord_update(p,p_prev,P,rand_coord,rand_coord_prev,kapa,lamdas,var,var,false);
}else
hit_and_run(p,P,var,var);

if (var.coordinate) {//Compute the first point for the CDHR
rand_coord = uidist(rng);
kapa = urdist(rng);
std::pair <FT, FT> bpair = P.line_intersect_coord(p, rand_coord, lamdas);
p_prev = p;
p.set_coord(rand_coord, p[rand_coord] + bpair.first + kapa * (bpair.second - bpair.first));
} else
hit_and_run(p, P, var, var);

for (int i = 1; i <= rnum; ++i) {

for (int j = 0; j < walk_len; ++j) {
if (var.coordinate) {
rand_coord_prev = rand_coord;
rand_coord = uidist(rng);
kapa = urdist(rng);
hit_and_run_coord_update(p, p_prev, P, rand_coord, rand_coord_prev, kapa, lamdas, var, var);
} else
hit_and_run(p, P, var, var);
}
randPoints.push_back(p);
//if(birk) birk_sym(P,randPoints,p);
}

//if(rand_only) std::cout<<p<<std::endl;
//if(print) std::cout<<"("<<i<<") Random point: "<<p<<std::endl;
}



template <class T, class K>
int rand_point_generator(BallIntersectPolytope<T> &PBLarge,
int rand_point_generator(BallIntersectPolytope<T,NT> &PBLarge,
Point &p, // a point to start
int rnum,
int walk_len,
K &randPoints,
BallIntersectPolytope<T> &PBSmall,
BallIntersectPolytope<T,NT> &PBSmall,
int &nump_PBSmall,
vars &var // constans for volume
)
{
//std::cout<<"EDW!: "<<std::endl;
vars &var) { // constans for volume
typedef typename Point::FT FT;
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::uniform_real_distribution<NT> urdist = var.urdist;
//std::uniform_int_distribution<int> uidist(0,n-1);

std::vector<NT> lamdas(PBLarge.num_of_hyperplanes(),NT(0));
//int rand_coord = rand()%n;
//double kapa = double(rand())/double(RAND_MAX);
int rand_coord = uidist(rng);
double kapa = urdist(rng);
boost::random::uniform_real_distribution<> urdist(0, 1);
boost::random::uniform_int_distribution<> uidist(0, n - 1);

std::vector <FT> lamdas(PBLarge.num_of_hyperplanes(), NT(0));
int rand_coord, rand_coord_prev;
FT kapa;
Point p_prev = p;

if(var.coordinate)
hit_and_run_coord_update(p,p_prev,PBLarge,rand_coord,rand_coord,kapa,lamdas,var,var,true);
else
hit_and_run(p,PBLarge,var,var);

for(int i=1; i<=rnum; ++i){
for(int j=0; j<walk_len; ++j){
int rand_coord_prev = rand_coord;
//rand_coord = rand()%n;
//kapa = double(rand())/double(RAND_MAX);
rand_coord = uidist(rng);
kapa = urdist(rng);
if(var.coordinate)
hit_and_run_coord_update(p,p_prev,PBLarge,rand_coord,rand_coord_prev,kapa,lamdas,var,var,false);
else
hit_and_run(p,PBLarge,var,var);
if (var.coordinate) {//Compute the first point for the CDHR
rand_coord = uidist(rng);
kapa = urdist(rng);
std::pair <FT, FT> bpair = PBLarge.line_intersect_coord(p, rand_coord, lamdas);
p_prev = p;
p.set_coord(rand_coord, p[rand_coord] + bpair.first + kapa * (bpair.second - bpair.first));
} else {
hit_and_run(p, PBLarge, var, var);
}

for (int i = 1; i <= rnum; ++i) {
for (int j = 0; j < walk_len; ++j) {
if (var.coordinate) {
rand_coord_prev = rand_coord;
rand_coord = uidist(rng);
kapa = urdist(rng);
hit_and_run_coord_update(p, p_prev, PBLarge, rand_coord, rand_coord_prev, kapa, lamdas, var, var);
} else
hit_and_run(p, PBLarge, var, var);
}
if(PBSmall.second().is_in(p) == -1){//is in
if (PBSmall.second().is_in(p) == -1) {//is in
randPoints.push_back(p);
++nump_PBSmall;
}
}
//if(rand_only) std::cout<<p<<std::endl;
//if(print) std::cout<<"("<<i<<") Random point: "<<p<<std::endl;
}

// ----- HIT AND RUN FUNCTIONS ------------ //
@@ -211,71 +197,42 @@ template <class T>
int hit_and_run(Point &p,
T &P,
vars &var,
vars &var2)
{
vars &var2) {
typedef typename Point::FT FT;
int n = var.n;
double err = var.err;
RNGType &rng = var.rng;
//std::uniform_real_distribution<NT> &urdist = var.urdist;
boost::random::uniform_real_distribution<> urdist(0,1);
//std::uniform_real_distribution<NT> &urdist1 = var.urdist1;
boost::random::uniform_real_distribution<> urdist(0, 1);

Point origin(n);

Random_points_on_sphere_d<Point> gen (n, 1.0);
Point l = gen.sample_point(rng);// - CGAL::Origin();
//Point l2=origin;
//Vector b1 = line_bisect(p,l,P,var,var2);
//Vector b2 = line_bisect(p,-l,P,var,var2);
//std::pair<Point,Point> ppair = P.line_intersect(p,l);
std::pair<NT,NT> dbpair = P.line_intersect(p,l);
NT min_plus = dbpair.first;
NT max_minus = dbpair.second;
Point b1 = (min_plus*l)+p;
Point b2 = (max_minus*l)+p;
//Point b1 = ppair.first;// - origin;
//Point b2 = ppair.second;// - origin;
//std::cout<<"b1="<<b1<<"b2="<<b2<<std::endl;

NT lambda = urdist(rng);
p = (lambda*b1);
p=((1-lambda)*b2) + p;
Random_points_on_sphere_d<Point> gen(n, 1.0);
Point l = gen.sample_point(rng);

std::pair <FT, FT> dbpair = P.line_intersect(p, l);
FT min_plus = dbpair.first;
FT max_minus = dbpair.second;
Point b1 = (min_plus * l) + p;
Point b2 = (max_minus * l) + p;
FT lambda = urdist(rng);
p = (lambda * b1);
p = ((1 - lambda) * b2) + p;
return 1;
}


//hit-and-run with orthogonal directions and update
template <class T>
template <class T, typename FT>
int hit_and_run_coord_update(Point &p,
Point &p_prev,
T &P,
int rand_coord,
int rand_coord_prev,
double kapa,
std::vector<NT> &lamdas,
FT kapa,
std::vector<FT> &lamdas,
vars &var,
vars &var2,
bool init)
{
//std::cout<<"[1]P dim: "<<P.dimension()<<std::endl;
std::pair<NT,NT> bpair;
// EXPERIMENTAL
//if(var.NN)
// bpair = P.query_dual(p,rand_coord);
//else
bpair = P.line_intersect_coord(p,p_prev,rand_coord,rand_coord_prev,lamdas,init);
//std::cout<<"[2]P dim: "<<P.dimension()<<std::endl;
//std::cout<<"original:"<<bpair.first<<" "<<bpair.second<<std::endl;
//std::cout<<"-----------"<<std::endl;
//TODO: only change one coordinate of *r* avoid addition + construction
//std::vector<NT> v(P.dimension(),NT(0));
//v[rand_coord] = bpair.first + kapa * (bpair.second - bpair.first);
//Point vp(P.dimension(),v.begin(),v.end());
vars &var2) {
std::pair <FT, FT> bpair = P.line_intersect_coord(p, p_prev, rand_coord, rand_coord_prev, lamdas);
p_prev = p;
//std::cout<<"v dim: "<<v.size()<<std::endl;
// std::cout<<"P dim: "<<P.dimension()<<std::endl;
//p = p + vp;
p.set_coord(rand_coord , p[rand_coord] + bpair.first + kapa * (bpair.second - bpair.first));
p.set_coord(rand_coord, p[rand_coord] + bpair.first + kapa * (bpair.second - bpair.first));
return 1;
}

29 changes: 13 additions & 16 deletions include/volume.h
Original file line number Diff line number Diff line change
@@ -14,7 +14,6 @@
#include <algorithm>
#include <math.h>
#include <chrono>
//#include <random.h>
#include "cartesian_geom/cartesian_kernel.h"
#include "random.hpp"
#include "random/uniform_int.hpp"
@@ -27,6 +26,8 @@


typedef double NT;
//typedef float NT;
//typedef long double NT;
typedef Cartesian<NT> Kernel;
typedef Kernel::Point Point;
typedef boost::mt19937 RNGType; ///< mersenne twister generator
@@ -39,14 +40,12 @@ struct vars{
int n,
int walk_steps,
int n_threads,
const double err,
double error,
const NT err,
NT error,
const int lw,
double up,
NT up,
const int L,
RNGType &rng,
//std::uniform_real_distribution<NT> urdist,
//std::uniform_real_distribution<NT> urdist1,
boost::random::uniform_real_distribution<>(urdist),
boost::random::uniform_real_distribution<> urdist1,
bool verbose,
@@ -65,14 +64,12 @@ struct vars{
int n;
int walk_steps;
int n_threads;
const double err;
double error;
const NT err;
NT error;
const int lw;
double up;
NT up;
const int L;
RNGType &rng;
//std::uniform_real_distribution<NT> urdist;
//std::uniform_real_distribution<NT> urdist1;
boost::random::uniform_real_distribution<>(urdist);
boost::random::uniform_real_distribution<> urdist1;
bool verbose;
@@ -129,9 +126,9 @@ template <class T>
NT volume(T &P,
vars &var, // constans for volume
vars &var2, // constants for optimization in case of MinkSums
std::pair<Point,double> CheBall) //Chebychev ball
std::pair<Point,NT> CheBall) //Chebychev ball
{
typedef BallIntersectPolytope<T> BallPoly;
typedef BallIntersectPolytope<T,NT> BallPoly;

bool round = var.round;
bool print = var.verbose;
@@ -140,7 +137,7 @@ NT volume(T &P,
int rnum = var.m;
int walk_len = var.walk_steps;
int n_threads = var.n_threads;
const double err = var.err;
const NT err = var.err;
RNGType &rng = var.rng;
// Rotation: only for test with skinny polytopes and rounding
//std::cout<<"Rotate="<<rotate(P)<<std::endl;
@@ -149,7 +146,7 @@ NT volume(T &P,
//0. Rounding of the polytope if round=true
Point c=CheBall.first;
NT radius=CheBall.second;
double round_value=1;
NT round_value=1;
if(round){
if(print) std::cout<<"\nRounding.."<<std::endl;
double tstart1 = (double)clock()/(double)CLOCKS_PER_SEC;
@@ -187,7 +184,7 @@ NT volume(T &P,

// 4. Construct the sequence of balls
// 4a. compute the radius of the largest ball
double current_dist, max_dist=NT(0);
NT current_dist, max_dist=NT(0);
for(std::list<Point>::iterator pit=randPoints.begin(); pit!=randPoints.end(); ++pit){
current_dist=(*pit-c).squared_length();
if(current_dist>max_dist){
14 changes: 7 additions & 7 deletions test/chebychev_test.cpp
Original file line number Diff line number Diff line change
@@ -17,20 +17,20 @@ long int factorial(int n)


template <typename FilePath>
void cheb_test(FilePath f, double expected, double tolerance=0.001)
void cheb_test(FilePath f, NT expected, NT tolerance=0.001)
{
std::ifstream inp;
std::vector<std::vector<double> > Pin;
std::vector<std::vector<NT> > Pin;
inp.open(f,std::ifstream::in);
read_pointset(inp,Pin);
int n = Pin[0][1]-1;
Polytope<double> P;
Polytope<NT> P;
P.init(Pin);

// Setup the parameters
double walk_len=10 + n/10;
int walk_len=10 + n/10;
int nexp=1, n_threads=1;
double e=1, err=0.0000000001;
NT e=1, err=0.0000000001;
int rnum = std::pow(e,-2) * 400 * n * std::log(n);
RNGType rng(std::time(0));
boost::normal_distribution<> rdist(0,1);
@@ -43,15 +43,15 @@ void cheb_test(FilePath f, double expected, double tolerance=0.001)
//Compute chebychev ball//
std::cout << "\n--- Testing Chebchev ball computation of " << f << std::endl;
double tstart1 = (double)clock()/(double)CLOCKS_PER_SEC;
std::pair<Point,double> CheBall = P.chebyshev_center();
std::pair<Point,NT> CheBall = P.chebyshev_center();
double tstop1 = (double)clock()/(double)CLOCKS_PER_SEC;

//std::cout<<"Chebychev center is: "<<std::endl;

std::cout<<"\nradius is: "<<CheBall.second<<std::endl;
std::cout << "Chebychev time = " << tstop1 - tstart1 << std::endl;

double error = std::abs(CheBall.second-expected)/expected;
NT error = std::abs(CheBall.second-expected)/expected;

CHECK(error < tolerance);

18 changes: 9 additions & 9 deletions test/rounding_test.cpp
Original file line number Diff line number Diff line change
@@ -17,20 +17,20 @@ long int factorial(int n)


template <typename FilePath>
void rounding_test(FilePath f, bool rot, double expected, double tolerance=0.1)//, double expected, double tolerance=0.1)
void rounding_test(FilePath f, bool rot, NT expected, NT tolerance=0.1)
{
std::ifstream inp;
std::vector<std::vector<double> > Pin;
std::vector<std::vector<NT> > Pin;
inp.open(f,std::ifstream::in);
read_pointset(inp,Pin);
int n = Pin[0][1]-1;
Polytope<double> P;
Polytope<NT> P;
P.init(Pin);

// Setup the parameters
double walk_len=10 + n/10;
int walk_len=10 + n/10;
int nexp=1, n_threads=1;
double e=1, err=0.0000000001;
NT e=1, err=0.0000000001;
int rnum = std::pow(e,-2) * 400 * n * std::log(n);
RNGType rng(std::time(0));
boost::normal_distribution<> rdist(0,1);
@@ -42,13 +42,13 @@ void rounding_test(FilePath f, bool rot, double expected, double tolerance=0.1)/

//Compute chebychev ball//
std::cout << "\n--- Testing rounding of " << f << std::endl;
double rot_val;
NT rot_val;
if(rot){
std::cout << "\n--- Rotation is ON "<< std::endl;
rot_val = rotating(P);
std::cout << "Rotation value = "<<rot_val<<std::endl;
}
std::pair<Point,double> CheBall;// = solveLP(P.get_matrix(), P.dimension());
std::pair<Point,NT> CheBall;// = solveLP(P.get_matrix(), P.dimension());
Point c;//=CheBall.first;
NT radius;//=CheBall.second;
NT round_value=1.0, ratio1,ratio2;
@@ -81,13 +81,13 @@ void rounding_test(FilePath f, bool rot, double expected, double tolerance=0.1)/
//c=CheBall.first;
//radius=CheBall.second;

double vol = 0;
NT vol = 0;
unsigned int const num_of_exp = 10;
for (unsigned int i=0; i<num_of_exp; i++)
{
vol += round_value*volume(P,var,var,CheBall);
}
double error = std::abs(((vol/num_of_exp)-expected))/expected;
NT error = std::abs(((vol/num_of_exp)-expected))/expected;
std::cout << "Computed volume (average) = " << vol/num_of_exp << std::endl;
std::cout << "Expected volume = " << expected << std::endl;
CHECK(error < tolerance);
34 changes: 17 additions & 17 deletions test/vol.cpp
Original file line number Diff line number Diff line change
@@ -41,8 +41,8 @@ int main(const int argc, const char** argv)
//Deafault values
int n, nexp=1, n_threads=1;
int walk_len;//to be defined after n
double e=1;
double exactvol(-1.0);
NT e=1;
NT exactvol(-1.0);
bool verbose=false,
rand_only=false,
round_only=false,
@@ -126,7 +126,7 @@ int main(const int argc, const char** argv)
file=true;
std::cout<<"Reading input from file..."<<std::endl;
std::ifstream inp;
std::vector<std::vector<double> > Pin;
std::vector<std::vector<NT> > Pin;
inp.open(argv[++i],std::ifstream::in);
read_pointset(inp,Pin);
//std::cout<<"d="<<Pin[0][1]<<std::endl;
@@ -168,7 +168,7 @@ int main(const int argc, const char** argv)

std::ifstream inp2;
inp2.open("order_polytope.ine",std::ifstream::in);
std::vector<std::vector<double> > Pin;
std::vector<std::vector<NT> > Pin;
read_pointset(inp2,Pin);

//std::cout<<"d="<<Pin[0][1]<<std::endl;
@@ -236,7 +236,7 @@ int main(const int argc, const char** argv)

//Compute chebychev ball//
double tstart1 = (double)clock()/(double)CLOCKS_PER_SEC;
std::pair<Point,double> CheBall = P.chebyshev_center();
std::pair<Point,NT> CheBall = P.chebyshev_center();
double tstop1 = (double)clock()/(double)CLOCKS_PER_SEC;
if(verbose) std::cout << "Chebychev time = " << tstop1 - tstart1 << std::endl;
if(verbose){
@@ -256,8 +256,8 @@ int main(const int argc, const char** argv)

/* CONSTANTS */
//error in hit-and-run bisection of P
const double err=0.0000000001;
const double err_opt=0.01;
const NT err=0.0000000001;
const NT err_opt=0.01;
//bounds for the cube
const int lw=0, up=10000, R=up-lw;

@@ -286,16 +286,16 @@ int main(const int argc, const char** argv)

//RUN EXPERIMENTS
int num_of_exp=nexp;
double sum=0, sum_time=0;
double min,max;
std::vector<double> vs;
double average, std_dev;
double sum_time=0;
NT min,max,sum=0;
std::vector<NT> vs;
NT average, std_dev;
double Chebtime, sum_Chebtime=double(0);
NT vol;

for(int i=0; i<num_of_exp; ++i){
std::cout<<"Experiment "<<i+1<<" ";
Polytope<double> P_to_test(P);
Polytope<NT> P_to_test(P);
tstart = (double)clock()/(double)CLOCKS_PER_SEC;

// Setup the parameters
@@ -304,20 +304,20 @@ int main(const int argc, const char** argv)

if(round_only){
// Round the polytope and exit
std::pair<double,double> res_round;
std::pair<NT,NT> res_round;
res_round = rounding_min_ellipsoid(P,CheBall,var);
double round_value = res_round.first;
NT round_value = res_round.first;
std::cout<<"\n--------------\nRounded polytope\nH-representation\nbegin\n"<<std::endl;
P.print();
std::cout<<"end\n--------------\n"<<std::endl;
}else{
// Estimate the volume
vol = volume(P_to_test,var,var,CheBall);
vol = volume(P,var,var,CheBall);
//if(rotate) vol = std::sqrt(vol);
//std::cout<<vol<<std::endl;
}

double v1 = vol;
NT v1 = vol;

tstop = (double)clock()/(double)CLOCKS_PER_SEC;
//double v2 = volume2(P,n,rnum,walk_len,err,rng,get_snd_rand,urdist,urdist1);
@@ -339,7 +339,7 @@ int main(const int argc, const char** argv)
//Compute Statistics
average=sum/(i+1);
std_dev=0;
for(std::vector<double>::iterator vit=vs.begin(); vit!=vs.end(); ++vit){
for(std::vector<NT>::iterator vit=vs.begin(); vit!=vs.end(); ++vit){
std_dev += std::pow(*vit - average,2);
}
std_dev = std::sqrt(std_dev/(i+1));
16 changes: 8 additions & 8 deletions test/volume_test.cpp
Original file line number Diff line number Diff line change
@@ -16,20 +16,20 @@ long int factorial(int n)
}

template <typename FilePath>
void test_volume(FilePath f, double expected, double tolerance=0.1)
void test_volume(FilePath f, NT expected, NT tolerance=0.1)
{
std::ifstream inp;
std::vector<std::vector<double> > Pin;
std::vector<std::vector<NT> > Pin;
inp.open(f,std::ifstream::in);
read_pointset(inp,Pin);
int n = Pin[0][1]-1;
Polytope<double> P;
Polytope<NT> P;
P.init(Pin);

// Setup the parameters
double walk_len=10 + n/10;
int walk_len=10 + n/10;
int nexp=1, n_threads=1;
double e=1, err=0.0000000001;
NT e=1, err=0.0000000001;
int rnum = std::pow(e,-2) * 400 * n * std::log(n);
RNGType rng(std::time(0));
boost::normal_distribution<> rdist(0,1);
@@ -40,17 +40,17 @@ void test_volume(FilePath f, double expected, double tolerance=0.1)
urdist,urdist1,false,false,false,false,false,true);

//Compute chebychev ball//
std::pair<Point,double> CheBall = P.chebyshev_center();
std::pair<Point,NT> CheBall = P.chebyshev_center();

// Estimate the volume
std::cout << "--- Testing volume of " << f << std::endl;
double vol = 0;
NT vol = 0;
unsigned int const num_of_exp = 10;
for (unsigned int i=0; i<num_of_exp; i++)
{
vol += volume(P,var,var,CheBall);
}
double error = std::abs(((vol/num_of_exp)-expected))/expected;
NT error = std::abs(((vol/num_of_exp)-expected))/expected;
std::cout << "Computed volume (average) = " << vol/num_of_exp << std::endl;
std::cout << "Expected volume = " << expected << std::endl;
CHECK(error < tolerance);

0 comments on commit 52da8e8

Please sign in to comment.