Skip to content

Commit

Permalink
Merge pull request #3 from iemiris/master
Browse files Browse the repository at this point in the history 
merging ioannis' version with few minor corrections
  • Loading branch information
@vissarion
vissarion authored May 30, 2017
2 parents 075929a + 45a95ef commit 96a66f5
Show file tree
Hide file tree
Showing 5 changed files with 100 additions and 92 deletions.
22 changes: 16 additions & 6 deletions README
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,10 @@
Approximate computation of volume of convex bodies.

By Vissarion Fisikopoulos, N.K. Univ. Athens, Greece, 2012-2016.
Main algorithm based on: I.Z. Emiris and V. Fisikopoulos, "Efficient random-walk methods for approximating polytope volume", In Proc. ACM Symposium on Computational Geometry 2014, pages 318-325.

You may redistribute or modify the software under the GNU Lesser General Public License as published by Free Software Foundation, either version 3 of the License, or (at your option) any later version. It is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY, see vol.cpp.


This is an open-source C++ software for computing an approximation of the volume of a polytope given as an intersection of halfspaces.

Expand All @@ -13,7 +20,6 @@ Follow the CGAL installation manual. It states that CGAL requires the Boost libr
$ cmake .
$ make


------------------
2. Compile sources
------------------
Expand Down Expand Up @@ -52,14 +58,14 @@ After successful compilation you can use the software by command line.

./vol -h

will display a help message about the program's available options.
will display a help message about the program's available options.

***Example***

To estimate the volume of the 10-dimensional hypercube first prepare the file cube.ine as follows:
To estimate the volume of the 10-dimensional hypercube first prepare the file cube10.ine as follows:

======================
cube_10.ine
cube10.ine
H-representation
begin
20 11 real
Expand Down Expand Up @@ -87,9 +93,13 @@ end
input_incidence
=======================

The run the following command:
Then run the following command:

./vol -f1 polytope_examples/cube10.ine

which returns 17 numbers:

d m #experiments exactvolOr-1 approx [.,.] #randPoints walkLength meanVol [minVol,maxVol] stdDev errorVsExact maxminDivergence time timeChebyshevBall

--------------------
5. Linear extensions
Expand Down Expand Up @@ -117,5 +127,5 @@ Then running

./vol -fle simple_poset.txt

will give an estimation of the number of linear extensions which is 3.
will give an estimation of the number of linear extensions, which is 3.

2 changes: 1 addition & 1 deletion examples/polytope_examples/cube10.ine
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
cube_10.ine
cube10.ine
H-representation
begin
20 11 real
Expand Down
113 changes: 56 additions & 57 deletions examples/vol.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,20 +36,20 @@ int factorial(int n)
// oracles.

int main(const int argc, const char** argv)
{
{
//Deafault values
int n, nexp=1, n_threads=1;
int n, nexp=1, n_threads=1;
int walk_len;//to be defined after n
double e=1;
double exactvol(-1.0);
bool verbose=false,
rand_only=false,
round_only=false,
file=false,
round=false,
NN=false,
user_walk_len=false,
linear_extensions=false,
double e=1;
double exactvol(-1.0);
bool verbose=false,
rand_only=false,
round_only=false,
file=false,
round=false,
NN=false,
user_walk_len=false,
linear_extensions=false,
birk=false,
rotate=false,
experiments=true,
Expand All @@ -63,27 +63,28 @@ int main(const int argc, const char** argv)
exit(-2);
}

//parse command line input vars
for(int i=1;i<argc;++i){
bool correct=false;
//parse command line input vars
for(int i=1;i<argc;++i){
bool correct=false;

if(!strcmp(argv[i],"-h")||!strcmp(argv[i],"--help")){
std::cerr<<
"Usage:\n"<<
"-v, --verbose \n"<<
"-rdhr : use random directions HnR, default is coordinate directions HnR\n"
"-rand, --rand_only : generates only random points\n"<<
"-f1, --file1 [filename type Ax<=b] [epsilon] [walk length] [threads] [num of experiments]\n"<<
//"-f2, --file2 [filename type Ax=b,x>=0] [epsilon] [walk length] [threads] [num of experiments]\n"<<
"-f1, --file1 [filename_type_Ax<=b] [epsilon] [walk_length] [threads] [num_of_experiments]\n"<<
//"-f2, --file2 [filename_type_Ax=b,x>=0] [epsilon] [walk_length] [threads] [num_of_experiments]\n"<<
"-fle, --filele : counting linear extensions of a poset\n"<<
//"-c, --cube [dimension] [epsilon] [walk length] [threads] [num of experiments]\n"<<
//"-c, --cube [dimension] [epsilon] [walk length] [threads] [num_of_experiments]\n"<<
"--exact : the exact volume\n"<<
"--cube : input polytope is a cube\n"<<
"-r, --round : enables rounding of the polytope as a preprocess\n"<<
"-ro, --round_only : does only rounding to the polytope\n"<<
"-e, --error [epsilon] : the goal error of approximation\n"<<
"-e, --error epsilon : the goal error of approximation\n"<<
"-w, --walk_len [walk_len] : the random walk length (default 10)\n"<<
"-exp [#exps] : number of experiments (default 1)\n"<<
"-t, --threads [#threads] : the number of threads to be used\n"<<
"-t, --threads #threads : the number of threads to be used\n"<<
"-ΝΝ : use Nearest Neighbor search to compute the boundary oracles\n"<<
"-birk_sym : use symmetry to compute more random points (only for Birkhoff polytopes)\n"<<
std::endl;
Expand All @@ -93,12 +94,12 @@ int main(const int argc, const char** argv)
exactvol = std::pow(2,n);
//exactvol = std::pow(2,n)/std::tgamma(n+1);//factorial of a natural number n is gamma(n+1)
correct=true;
}
}
if(!strcmp(argv[i],"--exact")){
exactvol = atof(argv[++i]);
correct=true;
}
if(!strcmp(argv[i],"-v")||!strcmp(argv[i],"--verbose")){
}
if(!strcmp(argv[i],"-v")||!strcmp(argv[i],"--verbose")){
verbose=true;
std::cout<<"Verbose mode\n";
correct=true;
Expand All @@ -124,7 +125,7 @@ int main(const int argc, const char** argv)
n = Pin[0][1]-1;
P.init(Pin);
if (verbose && P.num_of_hyperplanes()<100){
std::cout<<"Input polytope: "<<n<<std::endl;
std::cout<<"Input polytope: "<<n<<std::endl;
P.print();
}
correct=true;
Expand All @@ -142,7 +143,7 @@ int main(const int argc, const char** argv)
P.rref();
n=P.dimension();
//if (verbose && P.num_of_hyperplanes()<1000){
// std::cout<<"Input polytope: "<<n<<std::endl;
// std::cout<<"Input polytope: "<<n<<std::endl;
// P.print();
//}
correct=true;
Expand All @@ -166,7 +167,7 @@ int main(const int argc, const char** argv)
n = Pin[0][1]-1;
P.init(Pin);
//if (verbose && P.num_of_hyperplanes()<100){
std::cout<<"Input polytope: "<<n<<std::endl;
std::cout<<"Input polytope: "<<n<<std::endl;
//P.print();
//}
linear_extensions = true;
Expand Down Expand Up @@ -201,7 +202,7 @@ int main(const int argc, const char** argv)
P.dual(-1);
*/
std::cout<<"flann software is needed for this option. Experimental feature."
<<"Currently under development."<<std::endl;
<<"Currently under development."<<std::endl;
correct=true;
}
if(!strcmp(argv[i],"-ro")){
Expand All @@ -222,8 +223,7 @@ int main(const int argc, const char** argv)
"\', try "<<argv[0]<<" --help"<<std::endl;
exit(-2);
}

}
}//for i

// Set the number of random walk steps
if(!user_walk_len)
Expand Down Expand Up @@ -285,7 +285,6 @@ int main(const int argc, const char** argv)
vars var(rnum,n,walk_len,n_threads,err,0,0,0,0,rng,get_snd_rand,
urdist,urdist1,verbose,rand_only,round,NN,birk,coordinate);


if(round_only){
// Round the polytope and exit
double round_value = rounding(P,var,var);
Expand All @@ -306,38 +305,38 @@ int main(const int argc, const char** argv)
tstop = (double)clock()/(double)CLOCKS_PER_SEC;
//double v2 = volume2(P,n,rnum,walk_len,err,rng,get_snd_rand,urdist,urdist1);

// Statistics
sum+=v1;
if(i==0){max=v1;min=v1;}
if(v1>max) max=v1;
if(v1<min) min=v1;
vs.push_back(v1);
sum_time += tstop-tstart;
sum_Chebtime += Chebtime;
// Statistics
sum+=v1;
if(i==0){max=v1;min=v1;}
if(v1>max) max=v1;
if(v1<min) min=v1;
vs.push_back(v1);
sum_time += tstop-tstart;
sum_Chebtime += Chebtime;

//std::cout<<"\t vol= "<<v1<<"\t time= "<<tstop-tstart;
if(round)
//std::cout<<"\t vol= "<<v1<<"\t time= "<<tstop-tstart;
if(round)
std::cout<<" (rounding is ON)";
std::cout<<std::endl;
std::cout<<std::endl;

//Compute Statistics
average=sum/(i+1);
std_dev=0;
for(std::vector<double>::iterator vit=vs.begin(); vit!=vs.end(); ++vit){
//Compute Statistics
average=sum/(i+1);
std_dev=0;
for(std::vector<double>::iterator vit=vs.begin(); vit!=vs.end(); ++vit){
std_dev += std::pow(*vit - average,2);
}
std_dev = std::sqrt(std_dev/(i+1));
}
std_dev = std::sqrt(std_dev/(i+1));

std::cout.precision(7);
std::cout.precision(7);

//MEMORY USAGE
//struct proc_t usage;
//look_up_our_self(&usage);
//MEMORY USAGE
//struct proc_t usage;
//look_up_our_self(&usage);

//Print statistics
//std::cout<<"\nSTATISTICS:"<<std::endl;
if (!experiments){
std::cout
//Print statistics
//std::cout<<"\nSTATISTICS:"<<std::endl;
if (!experiments){
std::cout
<<"Dimension= "
<<n<<" "
//<<argv[]<<" "
Expand Down Expand Up @@ -379,8 +378,8 @@ int main(const int argc, const char** argv)
<<(exactvol-average)/exactvol<<" "
<<std::endl;
}
}else
std::cout
} else
std::cout
<<n<<" "
//<<argv[]<<" "
<<P.num_of_hyperplanes()<<" "
Expand Down
10 changes: 5 additions & 5 deletions include/polytopes.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -138,14 +138,14 @@ class stdHPolytope{

// print polytope in input format
int print() {
std::cout<<" "<<_A.size()<<" "<<_d+1<<" float"<<std::endl;
for(typename stdMatrix::iterator mit=_A.begin(); mit<_A.end(); ++mit){
for(typename stdCoeffs::iterator lit=mit->begin(); lit<mit->end() ; ++lit)
std::cout<<" "<<_A.size()<<" "<<_d+1<<" float"<<std::endl;
for(typename stdMatrix::iterator mit=_A.begin(); mit<_A.end(); ++mit){
for(typename stdCoeffs::iterator lit=mit->begin(); lit<mit->end() ; ++lit)
std::cout<<*lit<<" ";
std::cout<<std::endl;
}
return 0;
}
return 0;
}

/*
int is_in(stdCoeffs p) {
Expand Down
45 changes: 22 additions & 23 deletions include/vol_rand.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,6 +41,10 @@
#include <vector>
#include <iostream>

#ifndef BOOST_MATH_CONSTANTS_CONSTANTS_INCLUDED
#include <boost/math/constants/constants.hpp>
#endif // Ioannis Emiris

#include <Eigen/Eigen>
//#include <Eigen/Cholesky>

Expand Down Expand Up @@ -83,23 +87,23 @@ typedef boost::variate_generator< RNGType, boost::normal_distribution<> > gener
struct vars{
public:
vars( int m,
int n,
int walk_steps,
int n_threads,
const double err,
const double err_opt,
const int lw,
double up,
const int L,
RNGType &rng,
generator
&get_snd_rand,
boost::random::uniform_real_distribution<> urdist,
boost::random::uniform_real_distribution<> urdist1,
bool verbose,
bool rand_only,
bool round,
bool NN,
int n,
int walk_steps,
int n_threads,
const double err,
const double err_opt,
const int lw,
double up,
const int L,
RNGType &rng,
generator
&get_snd_rand,
boost::random::uniform_real_distribution<> urdist,
boost::random::uniform_real_distribution<> urdist1,
bool verbose,
bool rand_only,
bool round,
bool NN,
bool birk,
bool coordinate
) :
Expand Down Expand Up @@ -781,10 +785,6 @@ NT volume1_reuse_test(T &P,
}






/*************************************************
/* VOLUME with random COORDINATES hit and run
* Here we reuse the random points we generate
Expand Down Expand Up @@ -1247,8 +1247,7 @@ template <class T>
NT volume2(T &P,
vars &var)
{
typedef BallIntersectPolytope<T> BallPoly;

typedef BallIntersectPolytope<T> BallPoly;
int n = var.n;
int rnum = var.m;
int walk_len = var.walk_steps;
Expand Down

0 comments on commit 96a66f5

Please sign in to comment.