From 3ec8e416fe9477b43c7623de1602dc29238a12ff Mon Sep 17 00:00:00 2001 From: Vissarion Fisikopoulos Date: Sat, 4 Apr 2020 09:51:07 +0300 Subject: [PATCH] Update readme, add `doc` directory and tutorials (#70) --- R-proj/src/direct_sampling.cpp | 2 +- R-proj/src/sample_points.cpp | 19 +- README.md | 260 +----------------- doc/cpp_interface.md | 196 +++++++++++++ doc/r_interface.md | 47 ++++ include/annealing/gaussian_annealing.h | 7 +- include/annealing/hpoly_annealing.h | 2 +- include/annealing/ratio_estimation.h | 14 +- include/cartesian_geom/point.h | 195 +++++++------ include/convex_bodies/ball.h | 46 ++-- include/convex_bodies/ballintersectconvex.h | 16 +- include/convex_bodies/ellipsoids.h | 9 +- include/convex_bodies/hpolytope.h | 169 ++++++------ include/convex_bodies/vpolyintersectvpoly.h | 21 +- include/convex_bodies/vpolytope.h | 66 ++--- include/convex_bodies/zonoIntersecthpoly.h | 22 +- include/convex_bodies/zpolytope.h | 32 +-- include/generators/h_polytopes_gen.h | 7 +- .../generators/known_polytope_generators.h | 1 + include/lp_oracles/solve_lp.h | 19 +- include/lp_oracles/vpolyoracles.h | 1 + include/lp_oracles/zpolyoracles.h | 1 + include/rotating.h | 1 + include/rounding.h | 21 +- include/samplers/gaussian_samplers.h | 22 +- include/samplers/samplers.h | 92 ++++--- include/volume/cooling_balls.h | 7 +- include/volume/cooling_hpoly.h | 1 + include/volume/exact_vols.h | 2 + include/volume/volume.h | 11 +- test/CMakeLists.txt | 11 +- 31 files changed, 697 insertions(+), 623 deletions(-) create mode 100644 doc/cpp_interface.md create mode 100644 doc/r_interface.md diff --git a/R-proj/src/direct_sampling.cpp b/R-proj/src/direct_sampling.cpp index 195cb5467..2c947b14c 100644 --- a/R-proj/src/direct_sampling.cpp +++ b/R-proj/src/direct_sampling.cpp @@ -117,7 +117,7 @@ Rcpp::NumericMatrix direct_sampling(Rcpp::Nullable body = R_NilValue unsigned int jj = 0; for (typename std::list::iterator rpit = randPoints.begin(); rpit!=randPoints.end(); rpit++, jj++) { - RetMat.col(jj) = Eigen::Map(&(*rpit).get_coeffs()[0], (*rpit).dimension()); + RetMat.col(jj) = rpit->getCoefficients(); } return Rcpp::wrap(RetMat); } diff --git a/R-proj/src/sample_points.cpp b/R-proj/src/sample_points.cpp index 09e4d3bc3..4194bd709 100644 --- a/R-proj/src/sample_points.cpp +++ b/R-proj/src/sample_points.cpp @@ -197,6 +197,7 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P = R_NilValue } } + if (Rcpp::as(random_walk).containsElementNamed("nburns")) { nburns = Rcpp::as(Rcpp::as(random_walk)["nburns"]); if (nburns < 0) { @@ -233,7 +234,8 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P = R_NilValue HP.normalize(); if (gaussian) { StartingPoint = StartingPoint - mode; - HP.shift(Eigen::Map(&mode.get_coeffs()[0], mode.dimension())); + HP.shift(mode.getCoefficients()); + } break; } @@ -252,7 +254,7 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P = R_NilValue if (billiard && diam < 0.0) VP.comp_diam(diam); if (gaussian) { StartingPoint = StartingPoint - mode; - VP.shift(Eigen::Map(&mode.get_coeffs()[0], mode.dimension())); + VP.shift(mode.getCoefficients()); } break; } @@ -271,7 +273,7 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P = R_NilValue if (billiard && diam < 0.0) ZP.comp_diam(diam); if (gaussian) { StartingPoint = StartingPoint - mode; - ZP.shift(Eigen::Map(&mode.get_coeffs()[0], mode.dimension())); + ZP.shift(mode.getCoefficients()); } break; } @@ -296,7 +298,7 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P = R_NilValue } if (gaussian) { StartingPoint = StartingPoint - mode; - VPcVP.shift(Eigen::Map(&mode.get_coeffs()[0], mode.dimension())); + VPcVP.shift(mode.getCoefficients()); } break; } @@ -342,14 +344,17 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P = R_NilValue MT RetMat(dim, numpoints); unsigned int jj = 0; + for (typename std::list::iterator rpit = randPoints.begin(); rpit!=randPoints.end(); rpit++, jj++) { if (gaussian) { - RetMat.col(jj) = Eigen::Map(&(*rpit).get_coeffs()[0], (*rpit).dimension()) + - Eigen::Map(&mode.get_coeffs()[0], mode.dimension()); + + RetMat.col(jj) = rpit->getCoefficients() + mode.getCoefficients(); + } else { - RetMat.col(jj) = Eigen::Map(&(*rpit).get_coeffs()[0], (*rpit).dimension()); + RetMat.col(jj) = (*rpit).getCoefficients(); } } + return Rcpp::wrap(RetMat); } diff --git a/README.md b/README.md index 98dd0586b..2307d6253 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -**VolEsti** is a `C++` library for volume approximation and sampling of convex bodies (*e.g.* polytopes) with an `R` and limited `python` interface. +**VolEsti** is a `C++` library for volume approximation and sampling of convex bodies (*e.g.* polytopes) with an `R` and limited `python` interface. **VolEsti** is part of the [GeomScale](https://geomscale.github.io) project. [![CRAN status](https://www.r-pkg.org/badges/version/volesti)](https://cran.r-project.org/package=volesti) [![CRAN downloads](https://cranlogs.r-pkg.org/badges/volesti)](https://cran.r-project.org/package=volesti) @@ -11,259 +11,25 @@ [![CircleCI master](https://circleci.com/gh/GeomScale/volume_approximation/tree/master.svg?style=shield)](https://circleci.com/gh/GeomScale/volume_approximation/tree/master) [![CircleCI master](https://circleci.com/gh/GeomScale/volume_approximation/tree/develop.svg?style=shield)](https://circleci.com/gh/GeomScale/volume_approximation/tree/develop) -### R Interface ------------- +### Documentation -#### Install Rcpp package - -* Install package-dependencies: `Rcpp`, `RcppEigen`, `BH`, `lpSolveAPI`. - -1. Then use devtools package to install `volesti` Rcpp package. In folder /root/R-proj Run: -```r -Rcpp::compileAttributes() -library(devtools) -devtools::build() -devtools::install() -library(volesti) -``` -2. You can use Rstudio as well to open `volesti.Rproj` and then click `build source Package` and then `Install and Restart` in `Build` at the menu bar. - -#### Generate CRAN version - -To generate the CRAN version of the R package follow the instructions below: - -1. From the command line navigate to folder `/cran_gen`. Then Run: -```r -source('genCRANpkg.R') -``` - -2. Open genCRANpkg.R script with `Rstudio` and run it. - -#### Run volesti from `R` -* The main function is `volume()`. It can be used to approximate the volume of a convex polytope given as a set of linear inequalities or a set of vertices (d-dimensional points) or as a Minkowski sum of segments (zonotope). There are two algorithms that can be used. The first is `SequenceOfBalls` and the second is `CoolingGaussian` (see References). -* The function `sample_points()` can be used to sample points from a convex polytope with uniform or spherical gaussian target distribution. -* The function `round_polytope()` can be used to round a convex polytope. -* The function `rand_rotate()` can be used to apply a random rotation to a convex polytope. - -For more details you can read the documentation in folder `/inst/doc`. - -#### Create pdf documentation from Rd files -* Install volesti library. -* In `R` mode (or in Rstudio) Run -``` -pack = "volesti" -path = find.package(pack) -system(paste(shQuote(file.path(R.home("bin"), "R")), - "CMD", "Rd2pdf", shQuote(path))) -``` -* The pdf will be created and saved in R-proj folder. -* We give such a documentation in /R-proj/doc folder. - -### C++ Interface ------------- - -#### Compile C++ sources and run tests - -To compile the C++ code you need the [lp_solve](http://lpsolve.sourceforge.net/5.5/) library. For example, for Unix/Linux you need `liblpsolve55.so` (this is available from the library's [webpage](http://lpsolve.sourceforge.net/5.5/) as well as a package in several linux distributions e.g. [debian](https://packages.debian.org/stretch/liblpsolve55-dev)). You have to specify the path to `liblpsolve55.so`, by running, in folder test: -``` -cmake -DLP_SOLVE=_PATH_TO_LIB_FILE_ . -make -``` -For example: `-DLP_SOLVE=/usr/lib/lpsolve/liblpsolve55.so` - -You can run the tests by `cmake test` or `ctest -jK` where `K` the number of `CPU` threads. By adding the option `--verbose` to `ctest` you get more information about the tests, *e.g.* time per test, volume computed and the name of the polytope or convex body. - -#### Polytope input - -The current version of the software assumes that the polytope is given in the form of linear inequalities i.e. {x \in R^d : Ax <= b} where A is a matrix of dimension m *x* d and b a vector of dimension m or as a set of m vertices {\in R^d} or as a Minkowski sum of m segments {\in R^d}. The input is described in an `.ine`-file (H-polytopes) or in a `.ext` file (V-polytopes or zonotopes). The `.ine` file is described as follows: - -``` -various comments -begin -m d+1 numbertype -b -A -end -various options -``` - -The `.ext` file is described as follows: -``` -various comments -begin -m d+1 numbertype -1 v_1 -.. ... -1 v_m -end -various options -``` -In V-polytope case v_i are vertices and in zonotope case they are segments. - -This filestype (or similar) is used by a number of other software in polyhedral computation (e.g. `cdd`, `vinci`, `latte`). In the current version of the software, the options are treated as comments and the numbertype as C++ double type. -If your input has equality constraints then you have to transform it in the form that only contain linear inequalities which described above by using some other software. We recommend to use latte https://www.math.ucdavis.edu/~latte for this transformation. - -#### Run volesti from command line - -After successful compilation you can use the software by command line. For example, the following command `./vol -h` 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 `cube10.ine` as follows: - -``` -cube10.ine -H-representation -begin - 20 11 real - 1 1 0 0 0 0 0 0 0 0 0 - 1 0 1 0 0 0 0 0 0 0 0 - 1 0 0 1 0 0 0 0 0 0 0 - 1 0 0 0 1 0 0 0 0 0 0 - 1 0 0 0 0 1 0 0 0 0 0 - 1 0 0 0 0 0 1 0 0 0 0 - 1 0 0 0 0 0 0 1 0 0 0 - 1 0 0 0 0 0 0 0 1 0 0 - 1 0 0 0 0 0 0 0 0 1 0 - 1 0 0 0 0 0 0 0 0 0 1 - 1 -1 0 0 0 0 0 0 0 0 0 - 1 0 -1 0 0 0 0 0 0 0 0 - 1 0 0 -1 0 0 0 0 0 0 0 - 1 0 0 0 -1 0 0 0 0 0 0 - 1 0 0 0 0 -1 0 0 0 0 0 - 1 0 0 0 0 0 -1 0 0 0 0 - 1 0 0 0 0 0 0 -1 0 0 0 - 1 0 0 0 0 0 0 0 -1 0 0 - 1 0 0 0 0 0 0 0 0 -1 0 - 1 0 0 0 0 0 0 0 0 0 -1 -end -input_incidence -``` - -Then to use SequenceOfBalls (SOB) algorithm run the following command: -``` -./vol -f1 cube_10.ine -``` - -which returns 17 numbers: -```d m #experiments exactvolOr-1 approxVolume [.,.] #randPoints walkLength meanVol [minVol,maxVol] stdDev errorVsExact maxminDivergence time timeChebyshevBall``` - -To use CoolingGaussian (CG) algorithm run the following command: -``` -./vol -f1 cube_10.ine -cg -``` -which returns the same output as before. - -To estimate the volume of a 10-dimensional V-cross polytope described in `cross_10.ext` as follows: -``` -cross_10.ext -V-representation -begin - 20 11 integer - 1 1 0 0 0 0 0 0 0 0 0 - 1 0 1 0 0 0 0 0 0 0 0 - 1 0 0 1 0 0 0 0 0 0 0 - 1 0 0 0 1 0 0 0 0 0 0 - 1 0 0 0 0 1 0 0 0 0 0 - 1 0 0 0 0 0 1 0 0 0 0 - 1 0 0 0 0 0 0 1 0 0 0 - 1 0 0 0 0 0 0 0 1 0 0 - 1 0 0 0 0 0 0 0 0 1 0 - 1 0 0 0 0 0 0 0 0 0 1 - 1 -1 0 0 0 0 0 0 0 0 0 - 1 0 -1 0 0 0 0 0 0 0 0 - 1 0 0 -1 0 0 0 0 0 0 0 - 1 0 0 0 -1 0 0 0 0 0 0 - 1 0 0 0 0 -1 0 0 0 0 0 - 1 0 0 0 0 0 -1 0 0 0 0 - 1 0 0 0 0 0 0 -1 0 0 0 - 1 0 0 0 0 0 0 0 -1 0 0 - 1 0 0 0 0 0 0 0 0 -1 0 - 1 0 0 0 0 0 0 0 0 0 -1 -end -hull -incidence -``` -Run: -``` -./vol -f2 cross_10.ext -``` -which returns the same output as before. - -To estimate the volume of a 4-dimensional zonotope defined by the Minkowski sum of 8 segments described in `zonotope_4_8.ext` as follows: -``` -zonotope_4_8.ext -Zonotpe -begin - 8 5 real - 1 0.981851 -0.188734 -0.189761 0.0812645 - 1 -0.0181493 0.811266 -0.189761 0.0812645 - 1 -0.0181493 -0.188734 0.810239 0.0812645 - 1 -0.0181493 -0.188734 -0.189761 1.08126 - 1 -0.177863 0.437661 -0.0861379 -0.674634 - 1 0.737116 -0.204646 -0.540973 -0.471883 - 1 -0.684154 0.262324 0.292341 -0.265955 - 1 -0.802502 -0.740403 0.0938152 0.0874131 -end -hull -incidence -``` -Run: -``` -./vol -f3 zonotope_4_8.ext -``` -Flag `-v` enables the print mode. - -#### Generate polytopes - -You can use executable `generator` to generate polytopes (hypercubes, simplices, cross polytopes, skinny hypercubes (only in H-representation), product of two simplices (only in H-representation) and zonotoes. For example: - -1. To generate a 10-dimensional hypercube in H-representation run: -``` -./generate -cube -h -d 10 -``` - -2. To generate a 20-dimensional simplex in V-representaion run: -``` -./generate -simplex -v -d 20 -``` - -3. To generate a 5-dimensional zonotope defined by 10 segments run: -``` -./generate -zonotope -d 5 -m 10 -``` - -Command `./generate -help` will display a help message about the program's available options. - -#### Sampling - -You can sample from a convex polytope uniformly or from the spherical gaussian distribution. For example: - -1. To sample uniformly from the 10-dimensional hypercube, run: -``` -./vol -f1 cube_10.ine -rand -nsample 1000 -``` -Flag -nsample declares the number of points we wish to sample (default is 100). - -2. To sample from the gaussian distribution, run: -``` -./vol -f1 cube_10.ine -rand -nsample 1300 -gaussian -variance 1.5 -``` -Flag `-variance` declares the variance (default is 1.0). The center of the spherical gaussian is the Chebychev center for H-polytopes, or the origin for zonotopes. For V-polytopes is the chebychev center of the simplex that is defined by a random choice of d+1 vertices. - -3. To sample from a zonotope described in zonotope.ext file run: -``` -./vol -f3 zonotope.ext -rand -nsample 1500 -``` -For V-polytopes use flag `-f2` before the `.ext` file. In all cases use flag `-v` to print the excecutional time. +* [Using the R interface](doc/r_interface.md) +* [Using the C++ Interface](doc/cpp_interface.md) +* [Wikipage with tutorials and demos](https://github.com/GeomScale/volume_approximation/wiki) +* [Tutorial given to PyData meetup](https://vissarion.github.io/tutorials/volesti_tutorial_pydata.html) ### Credits -Copyright (c) 2012-2018 Vissarion Fisikopoulos -Copyright (c) 2018 Apostolos Chalkis +Copyright (c) 2012-2020 Vissarion Fisikopoulos +Copyright (c) 2018-2020 Apostolos Chalkis 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. -Main development by Vissarion Fisikopoulos while he was affiliated with University of Athens (UoA, Greece), University of Brussels (ULB, Belgium), and Chalkis Apostolos affiliated with University of Athens (UoA, Greece). +Main development started in 2012 by Vissarion Fisikopoulos while he was affiliated with University of Athens (UoA, Greece), and then with University of Brussels (ULB, Belgium). The main sampling and volume algorithms were developed at that time. + +Later, Chalkis Apostolos affiliated with University of Athens (UoA, Greece) and mentored by V.F. and partially supported by GSoC'2018, 2019 programs enhanced volesti with more algorithms for sampling and volume as well as an `R` interface. + +We acknowledge several contributions by the open-source community, most notably a `python` interface by Pedro Zuidberg Dos Martires affiliated with KU Leuven. ### Publications diff --git a/doc/cpp_interface.md b/doc/cpp_interface.md new file mode 100644 index 000000000..1296c98f3 --- /dev/null +++ b/doc/cpp_interface.md @@ -0,0 +1,196 @@ +### C++ Interface +------------ + +#### Compile C++ sources and run tests + +To compile the C++ code you need the [lp_solve](http://lpsolve.sourceforge.net/5.5/) library. For example, for Unix/Linux you need `liblpsolve55.so` (this is available from the library's [webpage](http://lpsolve.sourceforge.net/5.5/) as well as a package in several linux distributions e.g. [debian](https://packages.debian.org/stretch/liblpsolve55-dev)). You have to specify the path to `liblpsolve55.so`, by running, in folder test: +``` +cmake -DLP_SOLVE=_PATH_TO_LIB_FILE_ . +make +``` +For example: `-DLP_SOLVE=/usr/lib/lpsolve/liblpsolve55.so` + +You can run the tests by `cmake test` or `ctest -jK` where `K` the number of `CPU` threads. By adding the option `--verbose` to `ctest` you get more information about the tests, *e.g.* time per test, volume computed and the name of the polytope or convex body. + +#### Polytope input + +The current version of the software assumes that the polytope is given in the form of linear inequalities i.e. {x \in R^d : Ax <= b} where A is a matrix of dimension m *x* d and b a vector of dimension m or as a set of m vertices {\in R^d} or as a Minkowski sum of m segments {\in R^d}. The input is described in an `.ine`-file (H-polytopes) or in a `.ext` file (V-polytopes or zonotopes). The `.ine` file is described as follows: + +``` +various comments +begin +m d+1 numbertype +b -A +end +various options +``` + +The `.ext` file is described as follows: +``` +various comments +begin +m d+1 numbertype +1 v_1 +.. ... +1 v_m +end +various options +``` +In V-polytope case v_i are vertices and in zonotope case they are segments. + +This filestype (or similar) is used by a number of other software in polyhedral computation (e.g. `cdd`, `vinci`, `latte`). In the current version of the software, the options are treated as comments and the numbertype as C++ double type. +If your input has equality constraints then you have to transform it in the form that only contain linear inequalities which described above by using some other software. We recommend to use latte https://www.math.ucdavis.edu/~latte for this transformation. + +#### Run volesti from command line + +After successful compilation you can use the software by command line. For example, the following command `./vol -h` 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 `cube10.ine` as follows: + +``` +cube10.ine +H-representation +begin + 20 11 real + 1 1 0 0 0 0 0 0 0 0 0 + 1 0 1 0 0 0 0 0 0 0 0 + 1 0 0 1 0 0 0 0 0 0 0 + 1 0 0 0 1 0 0 0 0 0 0 + 1 0 0 0 0 1 0 0 0 0 0 + 1 0 0 0 0 0 1 0 0 0 0 + 1 0 0 0 0 0 0 1 0 0 0 + 1 0 0 0 0 0 0 0 1 0 0 + 1 0 0 0 0 0 0 0 0 1 0 + 1 0 0 0 0 0 0 0 0 0 1 + 1 -1 0 0 0 0 0 0 0 0 0 + 1 0 -1 0 0 0 0 0 0 0 0 + 1 0 0 -1 0 0 0 0 0 0 0 + 1 0 0 0 -1 0 0 0 0 0 0 + 1 0 0 0 0 -1 0 0 0 0 0 + 1 0 0 0 0 0 -1 0 0 0 0 + 1 0 0 0 0 0 0 -1 0 0 0 + 1 0 0 0 0 0 0 0 -1 0 0 + 1 0 0 0 0 0 0 0 0 -1 0 + 1 0 0 0 0 0 0 0 0 0 -1 +end +input_incidence +``` + +Then to use SequenceOfBalls (SOB) algorithm run the following command: +``` +./vol -f1 cube_10.ine +``` + +which returns 17 numbers: +```d m #experiments exactvolOr-1 approxVolume [.,.] #randPoints walkLength meanVol [minVol,maxVol] stdDev errorVsExact maxminDivergence time timeChebyshevBall``` + +To use CoolingGaussian (CG) algorithm run the following command: +``` +./vol -f1 cube_10.ine -cg +``` +which returns the same output as before. + +To estimate the volume of a 10-dimensional V-cross polytope described in `cross_10.ext` as follows: +``` +cross_10.ext +V-representation +begin + 20 11 integer + 1 1 0 0 0 0 0 0 0 0 0 + 1 0 1 0 0 0 0 0 0 0 0 + 1 0 0 1 0 0 0 0 0 0 0 + 1 0 0 0 1 0 0 0 0 0 0 + 1 0 0 0 0 1 0 0 0 0 0 + 1 0 0 0 0 0 1 0 0 0 0 + 1 0 0 0 0 0 0 1 0 0 0 + 1 0 0 0 0 0 0 0 1 0 0 + 1 0 0 0 0 0 0 0 0 1 0 + 1 0 0 0 0 0 0 0 0 0 1 + 1 -1 0 0 0 0 0 0 0 0 0 + 1 0 -1 0 0 0 0 0 0 0 0 + 1 0 0 -1 0 0 0 0 0 0 0 + 1 0 0 0 -1 0 0 0 0 0 0 + 1 0 0 0 0 -1 0 0 0 0 0 + 1 0 0 0 0 0 -1 0 0 0 0 + 1 0 0 0 0 0 0 -1 0 0 0 + 1 0 0 0 0 0 0 0 -1 0 0 + 1 0 0 0 0 0 0 0 0 -1 0 + 1 0 0 0 0 0 0 0 0 0 -1 +end +hull +incidence +``` +Run: +``` +./vol -f2 cross_10.ext +``` +which returns the same output as before. + +To estimate the volume of a 4-dimensional zonotope defined by the Minkowski sum of 8 segments described in `zonotope_4_8.ext` as follows: +``` +zonotope_4_8.ext +Zonotpe +begin + 8 5 real + 1 0.981851 -0.188734 -0.189761 0.0812645 + 1 -0.0181493 0.811266 -0.189761 0.0812645 + 1 -0.0181493 -0.188734 0.810239 0.0812645 + 1 -0.0181493 -0.188734 -0.189761 1.08126 + 1 -0.177863 0.437661 -0.0861379 -0.674634 + 1 0.737116 -0.204646 -0.540973 -0.471883 + 1 -0.684154 0.262324 0.292341 -0.265955 + 1 -0.802502 -0.740403 0.0938152 0.0874131 +end +hull +incidence +``` +Run: +``` +./vol -f3 zonotope_4_8.ext +``` +Flag `-v` enables the print mode. + +#### Generate polytopes + +You can use executable `generator` to generate polytopes (hypercubes, simplices, cross polytopes, skinny hypercubes (only in H-representation), product of two simplices (only in H-representation) and zonotoes. For example: + +1. To generate a 10-dimensional hypercube in H-representation run: +``` +./generate -cube -h -d 10 +``` + +2. To generate a 20-dimensional simplex in V-representaion run: +``` +./generate -simplex -v -d 20 +``` + +3. To generate a 5-dimensional zonotope defined by 10 segments run: +``` +./generate -zonotope -d 5 -m 10 +``` + +Command `./generate -help` will display a help message about the program's available options. + +#### Sampling + +You can sample from a convex polytope uniformly or from the spherical gaussian distribution. For example: + +1. To sample uniformly from the 10-dimensional hypercube, run: +``` +./vol -f1 cube_10.ine -rand -nsample 1000 +``` +Flag -nsample declares the number of points we wish to sample (default is 100). + +2. To sample from the gaussian distribution, run: +``` +./vol -f1 cube_10.ine -rand -nsample 1300 -gaussian -variance 1.5 +``` +Flag `-variance` declares the variance (default is 1.0). The center of the spherical gaussian is the Chebychev center for H-polytopes, or the origin for zonotopes. For V-polytopes is the chebychev center of the simplex that is defined by a random choice of d+1 vertices. + +3. To sample from a zonotope described in zonotope.ext file run: +``` +./vol -f3 zonotope.ext -rand -nsample 1500 +``` +For V-polytopes use flag `-f2` before the `.ext` file. In all cases use flag `-v` to print the excecutional time. diff --git a/doc/r_interface.md b/doc/r_interface.md new file mode 100644 index 000000000..11cf12b18 --- /dev/null +++ b/doc/r_interface.md @@ -0,0 +1,47 @@ +### R Interface +------------ + +#### Install Rcpp package + +* Install package-dependencies: `Rcpp`, `RcppEigen`, `BH`, `lpSolveAPI`. + +1. Then use devtools package to install `volesti` Rcpp package. In folder /root/R-proj Run: +```r +Rcpp::compileAttributes() +library(devtools) +devtools::build() +devtools::install() +library(volesti) +``` +2. You can use Rstudio as well to open `volesti.Rproj` and then click `build source Package` and then `Install and Restart` in `Build` at the menu bar. + +#### Generate CRAN version + +To generate the CRAN version of the R package follow the instructions below: + +1. From the command line navigate to folder `/cran_gen`. Then Run: +```r +source('genCRANpkg.R') +``` + +2. Open genCRANpkg.R script with `Rstudio` and run it. + +#### Run volesti from `R` +* The main function is `volume()`. It can be used to approximate the volume of a convex polytope given as a set of linear inequalities or a set of vertices (d-dimensional points) or as a Minkowski sum of segments (zonotope). There are two algorithms that can be used. The first is `SequenceOfBalls` and the second is `CoolingGaussian` (see References). +* The function `sample_points()` can be used to sample points from a convex polytope with uniform or spherical gaussian target distribution. +* The function `round_polytope()` can be used to round a convex polytope. +* The function `rand_rotate()` can be used to apply a random rotation to a convex polytope. + +For more details you can read the documentation in folder `/inst/doc`. + +#### Create pdf documentation from Rd files +* Install volesti library. +* In `R` mode (or in Rstudio) Run +``` +pack = "volesti" +path = find.package(pack) +system(paste(shQuote(file.path(R.home("bin"), "R")), + "CMD", "Rd2pdf", shQuote(path))) +``` +* The pdf will be created and saved in R-proj folder. +* We give such a documentation in /R-proj/doc folder. diff --git a/include/annealing/gaussian_annealing.h b/include/annealing/gaussian_annealing.h index 01a1a0d33..eea190670 100644 --- a/include/annealing/gaussian_annealing.h +++ b/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 lamdas(P.num_of_hyperplanes(), NT(0)); + typedef Eigen::Matrix 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){ diff --git a/include/annealing/hpoly_annealing.h b/include/annealing/hpoly_annealing.h index 0a9f929bb..d262d84ee 100644 --- a/include/annealing/hpoly_annealing.h +++ b/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::iterator rpit=randPoints.begin(); NT max_norm = 0.0, iter_norm; for ( ; rpit!=randPoints.end(); rpit++) { - iter_norm = (G*Eigen::Map(&(*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); diff --git a/include/annealing/ratio_estimation.h b/include/annealing/ratio_estimation.h index efe5577c0..d8256fb60 100644 --- a/include/annealing/ratio_estimation.h +++ b/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::lowest(), max_val = std::numeric_limits::max(), val, lambda; size_t totCount = Ntot, countIn = Ntot * ratio; - std::vector last_W(W), lamdas(Pb1.num_of_hyperplanes()), Av(Pb1.num_of_hyperplanes()); + std::vector last_W(W); + typedef Eigen::Matrix VT; + VT lamdas, Av; + lamdas.setZero(Pb1.num_of_hyperplanes()); + Av.setZero(Pb1.num_of_hyperplanes()); + std::list randPoints; typename std::vector::iterator minmaxIt; typename std::list::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 last_W(W), lamdas(Pb1.num_of_hyperplanes()), Av(Pb1.num_of_hyperplanes()); + std::vector last_W(W); + typedef Eigen::Matrix 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 = "< e * std::abs(*mit) || - std::abs(*mit - *pit) > e * std::abs(*pit)) + for (i=0 ; icoeffs(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(const point& p) const { + return coeffs.dot(p.getCoefficients()); + } - FT dot(point& p){ - FT res=FT(0.0); - - typename Coeff::iterator pit=p.iter_begin(); - typename Coeff::iterator mit=coeffs.begin(); - for( ; pitcoeffs.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 -point operator* (typename K::FT const& k, point& p) { +point operator* (const typename K::FT& k, point& p) { return p * k; } #endif + + diff --git a/include/convex_bodies/ball.h b/include/convex_bodies/ball.h index 1a470e6bb..91a7cfb82 100644 --- a/include/convex_bodies/ball.h +++ b/include/convex_bodies/ball.h @@ -3,6 +3,8 @@ // Copyright (c) 20012-2020 Vissarion Fisikopoulos // Copyright (c) 2018-2020 Apostolos Chalkis +//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. + // Licensed under GNU LGPL.3, see LICENCE file #ifndef BALL_H @@ -15,6 +17,7 @@ struct Ball{ typedef Point BallPoint; typedef typename Point::FT NT; typedef typename std::vector::iterator viterator; + typedef Eigen::Matrix VT; Ball() {} @@ -44,26 +47,22 @@ struct Ball{ std::pair line_intersect(Point &r, Point &v) { - viterator vit=v.iter_begin(), cit=c.iter_begin(), rcit=r.iter_begin(); NT vrc(0), v2(0), rc2(0); - for( ; cit < c.iter_end() ; ++cit, ++rcit, ++vit){ - vrc += *vit * (*rcit); - v2 += *vit * (*vit); - rc2 += *rcit * (*rcit); - } + + vrc = v.dot(r); + v2 = v.dot(v); + rc2 = r.dot(r); NT disc_sqrt = std::sqrt(std::pow(vrc,2) - v2 * (rc2 - R)); return std::pair ((NT(-1)*vrc + disc_sqrt)/v2, (NT(-1)*vrc - disc_sqrt)/v2); } - std::pair line_intersect(Point &r, Point &v, const std::vector &Ar, - const std::vector &Av){ + std::pair line_intersect(Point &r, Point &v, const VT &Ar, const VT &Av){ return line_intersect(r, v); } - std::pair line_intersect(Point &r, Point &v, const std::vector &Ar, - const std::vector &Av, NT &lambda_prev) { + std::pair line_intersect(Point &r, Point &v, const VT &Ar, const VT &Av, NT &lambda_prev) { return line_intersect(r, v); } @@ -71,24 +70,22 @@ struct Ball{ return std::pair(line_intersect(r, v).first, 0); } - std::pair line_positive_intersect(Point &r, Point &v, const std::vector &Ar, - const std::vector &Av){ + std::pair line_positive_intersect(Point &r, Point &v, const VT &Ar, + const VT &Av){ return line_positive_intersect(r, v); } - std::pair line_positive_intersect(Point &r, Point &v, const std::vector &Ar, - const std::vector &Av, NT &lambda_prev){ + std::pair line_positive_intersect(Point &r, Point &v, const VT &Ar, + const VT &Av, NT &lambda_prev){ return line_positive_intersect(r, v); } std::pair line_intersect_coord(Point &r, const unsigned int &rand_coord) { - viterator rcit=r.iter_begin(); - NT vrc = *(rcit + rand_coord); + NT vrc = r[rand_coord]; NT rc2(R); - for( ; rcit < r.iter_end() ; ++rcit){ - rc2 -= *rcit * (*rcit); - } + rc2 -= r.dot(r); + NT disc_sqrt = std::sqrt(std::pow(vrc,2) + rc2); return std::pair (NT(-1)*vrc + disc_sqrt, NT(-1)*vrc - disc_sqrt); @@ -96,7 +93,7 @@ struct Ball{ } std::pair line_intersect_coord(Point &r, const unsigned int &rand_coord, - const std::vector &lamdas) { + const VT &lamdas) { return line_intersect_coord(r, rand_coord); } @@ -104,7 +101,7 @@ struct Ball{ const Point &r_prev, const unsigned int &rand_coord, const unsigned int &rand_coord_prev, - const std::vector &lamdas) { + const VT &lamdas) { return line_intersect_coord(r, rand_coord); } @@ -115,10 +112,9 @@ struct Ball{ void compute_reflection (Point &v, const Point &p, const int &facet) { Point s = p; - s = s * (1.0 / std::sqrt(s.squared_length())); - s = ((-2.0 * v.dot(s)) * s); - v = s + v; - + s *= (1.0 / std::sqrt(s.squared_length())); + s *= (-2.0 * v.dot(s)); + v += s; } private: diff --git a/include/convex_bodies/ballintersectconvex.h b/include/convex_bodies/ballintersectconvex.h index cb1c52e61..e670c65ad 100644 --- a/include/convex_bodies/ballintersectconvex.h +++ b/include/convex_bodies/ballintersectconvex.h @@ -4,6 +4,7 @@ // Copyright (c) 2018-2020 Apostolos Chalkis //Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 program. +//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. // Licensed under GNU LGPL.3, see LICENCE file @@ -19,6 +20,7 @@ class BallIntersectPolytope { public: typedef typename CBall::NT NT; typedef typename CBall::BallPoint Point; + typedef Eigen::Matrix VT; BallIntersectPolytope() {} @@ -67,15 +69,15 @@ class BallIntersectPolytope { std::max(polypair.second, ballpair.second)); } - std::pair line_intersect(Point &r, Point &v, std::vector &Ar, - std::vector &Av) { + std::pair line_intersect(Point &r, Point &v, VT &Ar, + VT &Av) { std::pair polypair = P.line_intersect(r, v, Ar, Av); std::pair ballpair = B.line_intersect(r, v); return std::pair(std::min(polypair.first, ballpair.first), std::max(polypair.second, ballpair.second)); } - std::pair line_intersect(Point &r, Point &v, std::vector &Ar, std::vector &Av, NT &lambda_prev) { + std::pair line_intersect(Point &r, Point &v, VT &Ar, VT &Av, NT &lambda_prev) { std::pair polypair = P.line_intersect(r, v, Ar, Av, lambda_prev); std::pair ballpair = B.line_intersect(r, v); @@ -83,7 +85,7 @@ class BallIntersectPolytope { std::max(polypair.second, ballpair.second)); } - std::pair line_positive_intersect(Point &r, Point &v, std::vector &Ar, std::vector &Av) { + std::pair line_positive_intersect(Point &r, Point &v, VT &Ar, VT &Av) { std::pair polypair = P.line_positive_intersect(r, v, Ar, Av); std::pair ball_lambda = B.line_positive_intersect(r, v); @@ -95,7 +97,7 @@ class BallIntersectPolytope { } - std::pair line_positive_intersect(Point &r, Point &v, std::vector &Ar, std::vector &Av, + std::pair line_positive_intersect(Point &r, Point &v, VT &Ar, VT &Av, NT &lambda_prev) { std::pair polypair = P.line_positive_intersect(r, v, Ar, Av, lambda_prev); @@ -110,7 +112,7 @@ class BallIntersectPolytope { //First coordinate ray shooting intersecting convex body std::pair line_intersect_coord(Point &r, const unsigned int &rand_coord, - std::vector &lamdas) { + VT &lamdas) { std::pair polypair = P.line_intersect_coord(r, rand_coord, lamdas); std::pair ballpair = B.line_intersect_coord(r, rand_coord); @@ -123,7 +125,7 @@ class BallIntersectPolytope { const Point &r_prev, const unsigned int &rand_coord, const unsigned int &rand_coord_prev, - std::vector &lamdas) { + VT &lamdas) { std::pair polypair = P.line_intersect_coord(r, r_prev, rand_coord, rand_coord_prev, lamdas); std::pair ballpair = B.line_intersect_coord(r, rand_coord); diff --git a/include/convex_bodies/ellipsoids.h b/include/convex_bodies/ellipsoids.h index 21926003b..72ef9cc88 100644 --- a/include/convex_bodies/ellipsoids.h +++ b/include/convex_bodies/ellipsoids.h @@ -4,6 +4,7 @@ // Copyright (c) 2018 Apostolos Chalkis //Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 program. +//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. // Licensed under GNU LGPL.3, see LICENCE file @@ -38,13 +39,7 @@ class copula_ellipsoid{ } NT mat_mult(Point p) { - VT q(dim); - unsigned int i = 0; - viterator pit = p.iter_begin(); - for ( ; pit!=p.iter_end(); ++pit, ++i){ - q(i)=(*pit); - } - return q.transpose()*G*q; + return p.getCoefficients().transpose()*G*p.getCoefficients(); } }; diff --git a/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h index 50946ffa3..f172fe534 100644 --- a/include/convex_bodies/hpolytope.h +++ b/include/convex_bodies/hpolytope.h @@ -4,6 +4,7 @@ // Copyright (c) 2018 Apostolos Chalkis //Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 program. +//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. // Licensed under GNU LGPL.3, see LICENCE file @@ -250,12 +251,14 @@ class HPolytope{ int is_in(const Point &p) const { NT sum; int m = A.rows(); - for (int i = 0; i < m; i++) { - sum = b(i); - for (unsigned int j = 0; j < _d; j++) sum -= A(i, j) * p[j]; + const NT* b_data = b.data(); + for (int i = 0; i < m; i++) { //Check if corresponding hyperplane is violated - if (sum < NT(0)) return 0; + if (*b_data - A.row(i) * p.getCoefficients() < NT(0)) + return 0; + + b_data++; } return -1; } @@ -266,30 +269,31 @@ class HPolytope{ std::pair line_intersect(Point &r, Point &v) { NT lamda = 0, min_plus = NT(maxNT), max_minus = NT(minNT); - NT sum_nom, sum_denom; + VT sum_nom, sum_denom; //unsigned int i, j; unsigned int j; int m = num_of_hyperplanes(); - viterator rit, vit; + + + sum_nom.noalias() = b - A * r.getCoefficients(); + sum_denom.noalias() = A * v.getCoefficients(); + + NT* sum_nom_data = sum_nom.data(); + NT* sum_denom_data = sum_denom.data(); for (int i = 0; i < m; i++) { - sum_nom = b(i); - sum_denom = NT(0); - j = 0; - rit = r.iter_begin(); - vit = v.iter_begin(); - for ( ; rit != r.iter_end(); rit++, vit++, j++){ - sum_nom -= A(i, j) * (*rit); - sum_denom += A(i, j) * (*vit); - } - if (sum_denom == NT(0)) { + + if (*sum_denom_data == NT(0)) { //std::cout<<"div0"< 0) min_plus = lamda; if (lamda > max_minus && lamda < 0) max_minus = lamda; } + + sum_nom_data++; + sum_denom_data++; } return std::pair(min_plus, max_minus); } @@ -304,35 +308,34 @@ class HPolytope{ bool pos = false) { NT lamda = 0, min_plus = NT(maxNT), max_minus = NT(minNT); - NT sum_nom, sum_denom, mult; + VT sum_nom; + NT mult; //unsigned int i, j; unsigned int j; int m = num_of_hyperplanes(), facet; - viterator rit, vit, Ariter = Ar.begin(), Aviter = Av.begin(); - - for (int i = 0; i < m; i++, ++Ariter, ++Aviter) { - sum_nom = NT(0); - sum_denom = NT(0); - j = 0; - rit = r.iter_begin(); - vit = v.iter_begin(); - for ( ; rit != r.iter_end(); rit++, vit++, j++){ - sum_nom -= A(i, j) * (*rit); - sum_denom += A(i, j) * (*vit); - } - (*Ariter) = -sum_nom; - (*Aviter) = sum_denom; - sum_nom += b(i); - if (sum_denom == NT(0)) { + + Ar.noalias() = A * r.getCoefficients(); + sum_nom = b - Ar; + Av.noalias() = A * v.getCoefficients();; + + + NT* Av_data = Av.data(); + NT* sum_nom_data = sum_nom.data(); + + for (int i = 0; i < m; i++) { + if (*Av_data == NT(0)) { //std::cout<<"div0"< 0) { min_plus = lamda; if (pos) facet = i; }else if (lamda > max_minus && lamda < 0) max_minus = lamda; } + + Av_data++; + sum_nom_data++; } if (pos) return std::pair(min_plus, facet); return std::pair(min_plus, max_minus); @@ -346,31 +349,32 @@ class HPolytope{ bool pos = false) { NT lamda = 0, min_plus = NT(maxNT), max_minus = NT(minNT); - NT sum_nom, sum_denom, mult; + VT sum_nom; + NT mult; //unsigned int i, j; unsigned int j; int m = num_of_hyperplanes(), facet; - viterator vit, Ariter = Ar.begin(), Aviter = Av.begin(); - - for (int i = 0; i < m; i++, ++Ariter, ++Aviter) { - (*Ariter) += lambda_prev * (*Aviter); - sum_nom = b(i) - (*Ariter); - sum_denom = NT(0); - j = 0; - vit = v.iter_begin(); - for ( ; vit != v.iter_end(); vit++, j++) sum_denom += A(i, j) * (*vit); - - (*Aviter) = sum_denom; - if (sum_denom == NT(0)) { + + Ar.noalias() += lambda_prev*Av; + sum_nom = b - Ar; + Av.noalias() = A * v.getCoefficients(); + + NT* sum_nom_data = sum_nom.data(); + NT* Av_data = Av.data(); + + for (int i = 0; i < m; i++) { + if (*Av_data == NT(0)) { //std::cout<<"div0"< 0) { min_plus = lamda; if (pos) facet = i; }else if (lamda > max_minus && lamda < 0) max_minus = lamda; } + Av_data++; + sum_nom_data++; } if (pos) return std::pair(min_plus, facet); return std::pair(min_plus, max_minus); @@ -379,48 +383,48 @@ class HPolytope{ // compute intersection point of a ray starting from r and pointing to v // with polytope discribed by A and b - std::pair line_positive_intersect(Point &r, Point &v, std::vector &Ar, - std::vector &Av) { + std::pair line_positive_intersect(Point &r, Point &v, VT& Ar, + VT& Av) { return line_intersect(r, v, Ar, Av, true); } // compute intersection point of a ray starting from r and pointing to v // with polytope discribed by A and b - std::pair line_positive_intersect(Point &r, Point &v, std::vector &Ar, - std::vector &Av, const NT &lambda_prev) { + std::pair line_positive_intersect(Point &r, Point &v, VT& Ar, + VT& Av, const NT &lambda_prev) { return line_intersect(r, v, Ar, Av, lambda_prev, true); } //First coordinate ray intersecting convex polytope std::pair line_intersect_coord(Point &r, const unsigned int &rand_coord, - std::vector &lamdas) { + VT& lamdas) { NT lamda = 0, min_plus = NT(maxNT), max_minus = NT(minNT); - NT sum_nom, sum_denom; + VT sum_denom; unsigned int j; int m = num_of_hyperplanes(); - viterator rit; + + sum_denom = A.col(rand_coord); + lamdas.noalias() = b - A * r.getCoefficients(); + + NT* lamda_data = lamdas.data(); + NT* sum_denom_data = sum_denom.data(); for (int i = 0; i < m; i++) { - sum_nom = b(i); - sum_denom = A(i, rand_coord); - rit = r.iter_begin(); - j = 0; - for (; rit != r.iter_end(); rit++, j++) { - sum_nom -= A(i, j) * (*rit); - } - lamdas[i] = sum_nom; - if (sum_denom == NT(0)) { + + if (*sum_denom_data == NT(0)) { //std::cout<<"div0"< 0) min_plus = lamda; if (lamda > max_minus && lamda < 0) max_minus = lamda; } + lamda_data++; + sum_denom_data++; } return std::pair(min_plus, max_minus); } @@ -431,29 +435,28 @@ class HPolytope{ const Point &r_prev, const unsigned int rand_coord, const unsigned int rand_coord_prev, - std::vector &lamdas) { - - viterator lamdait = lamdas.begin(); + VT& lamdas) { + ; NT lamda = 0, min_plus = NT(maxNT), max_minus = NT(minNT); - NT sum_nom, sum_denom, c_rand_coord, c_rand_coord_prev; + int m = num_of_hyperplanes(); + lamdas.noalias() += A.col(rand_coord_prev)* (r_prev[rand_coord_prev] - r[rand_coord_prev]); + NT* data = lamdas.data(); + for (int i = 0; i < m; i++) { - sum_denom = b(i); - c_rand_coord = A(i, rand_coord); - c_rand_coord_prev = A(i, rand_coord_prev); + NT a = A(i, rand_coord); - *lamdait = *lamdait + c_rand_coord_prev * (r_prev[rand_coord_prev] - r[rand_coord_prev]); - if (c_rand_coord == NT(0)) { + if (a == NT(0)) { //std::cout<<"div0"< 0) min_plus = lamda; if (lamda > max_minus && lamda < 0) max_minus = lamda; } - ++lamdait; + data++; } return std::pair(min_plus, max_minus); } @@ -467,7 +470,7 @@ class HPolytope{ // shift polytope by a point c void shift(const VT &c){ - b = b - A*c; + b -= A*c; } @@ -506,11 +509,11 @@ class HPolytope{ void compute_reflection(Point &v, const Point &p, const int facet) { - VT a = A.row(facet); - Point s(_d, std::vector(&a[0], a.data()+a.cols()*a.rows())); - s = ((-2.0 * v.dot(s)) * s); - v = s + v; - +// VT a = A.row(facet); +// Point s(_d, std::vector(&a[0], a.data()+a.cols()*a.rows())); +// Point s(A.row(facet)); +// s *= (-2.0 * v.dot(s)); + v += -2 * v.dot(A.row(facet)) * A.row(facet); } void free_them_all() {} diff --git a/include/convex_bodies/vpolyintersectvpoly.h b/include/convex_bodies/vpolyintersectvpoly.h index e49c0d45d..c1bb9a3eb 100644 --- a/include/convex_bodies/vpolyintersectvpoly.h +++ b/include/convex_bodies/vpolyintersectvpoly.h @@ -4,6 +4,7 @@ // Copyright (c) 2018 Apostolos Chalkis //Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 program. +//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. // Licensed under GNU LGPL.3, see LICENCE file @@ -210,16 +211,16 @@ class IntersectionOfVpoly { // compute intersection point of ray starting from r and pointing to v // with the V-polytope - std::pair line_intersect(const Point &r, const Point &v, const std::vector &Ar, - const std::vector &Av) { + std::pair line_intersect(const Point &r, const Point &v, const VT &Ar, + const VT &Av) { return line_intersect(r, v); } // compute intersection point of ray starting from r and pointing to v // with the V-polytope - std::pair line_intersect(const Point &r, const Point &v, const std::vector &Ar, - const std::vector &Av, const NT &lambda) { + std::pair line_intersect(const Point &r, const Point &v, const VT &Ar, + const VT &Av, const NT &lambda) { return line_intersect(r, v); } @@ -234,14 +235,14 @@ class IntersectionOfVpoly { return std::pair(P2pair.first, 2); } - std::pair line_positive_intersect(const Point &r, const Point &v, const std::vector &Ar, - const std::vector &Av) { + std::pair line_positive_intersect(const Point &r, const Point &v, const VT &Ar, + const VT &Av) { return line_positive_intersect(r, v); } - std::pair line_positive_intersect(const Point &r, const Point &v, const std::vector &Ar, - const std::vector &Av, const NT &lambda_prev) { + std::pair line_positive_intersect(const Point &r, const Point &v, const VT &Ar, + const VT &Av, const NT &lambda_prev) { return line_positive_intersect(r, v);//, Ar, Av); } @@ -250,7 +251,7 @@ class IntersectionOfVpoly { // with the V-polytope std::pair line_intersect_coord(const Point &r, const unsigned int &rand_coord, - const std::vector &lamdas) { + const VT &lamdas) { std::pair P1pair = P1.line_intersect_coord(r, rand_coord, lamdas); std::pair P2pair = P2.line_intersect_coord(r, rand_coord, lamdas); return std::pair(std::min(P1pair.first, P2pair.first), @@ -264,7 +265,7 @@ class IntersectionOfVpoly { const Point &r_prev, const unsigned int &rand_coord, const unsigned int &rand_coord_prev, - const std::vector &lamdas) { + const VT &lamdas) { return line_intersect_coord(r, rand_coord, lamdas); } diff --git a/include/convex_bodies/vpolytope.h b/include/convex_bodies/vpolytope.h index cce85a7fd..14d0135f8 100644 --- a/include/convex_bodies/vpolytope.h +++ b/include/convex_bodies/vpolytope.h @@ -4,6 +4,7 @@ // Copyright (c) 2018 Apostolos Chalkis //Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 program. +//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. // Licensed under GNU LGPL.3, see LICENCE file @@ -165,15 +166,11 @@ class VPolytope{ } Point get_mean_of_vertices() { - std::vector vec(_d); - Point xc(_d), temp(_d); + Point xc(_d); for (int i = 0; i < num_of_vertices(); ++i) { - for (int j = 0; j < _d; ++j) vec[j] = V(i,j); - - temp = Point(_d, vec.begin(), vec.end()); - xc = xc + temp; + xc.add(V.row(i)); } - xc = xc * (1.0/NT(num_of_vertices())); + xc *= (1.0/NT(num_of_vertices())); return xc; } @@ -222,17 +219,15 @@ class VPolytope{ std::pair result; for (j=1; jgetCoefficients() - p0.getCoefficients(); } Bg = B; B = B.inverse(); for (i=0; i Ap(_d,randPoints.size()); typename std::list::iterator rpit=randPoints.begin(); - typename std::vector::iterator qit; + unsigned int i, j = 0; for ( ; rpit!=randPoints.end(); rpit++, j++) { - qit = (*rpit).iter_begin(); i=0; - for ( ; qit!=(*rpit).iter_end(); qit++, i++){ - Ap(i,j)=double(*qit); + const NT* point_data = rpit->getCoefficients().data(); + + for ( i=0; i < rpit->dimension(); i++){ + Ap(i,j)=double(*point_data); + point_data++; } } boost::numeric::ublas::matrix Q(_d, _d); @@ -355,16 +352,16 @@ class VPolytope{ // compute intersection point of ray starting from r and pointing to v // with the V-polytope - std::pair line_intersect(const Point &r, const Point &v, const std::vector &Ar, - const std::vector &Av) { + std::pair line_intersect(const Point &r, const Point &v, const VT &Ar, + const VT &Av) { return intersect_double_line_Vpoly(V, r, v, row, colno); } // compute intersection point of ray starting from r and pointing to v // with the V-polytope - std::pair line_intersect(const Point &r, const Point &v, const std::vector &Ar, - const std::vector &Av, const NT &lambda_prev) { + std::pair line_intersect(const Point &r, const Point &v, const VT &Ar, + const VT &Av, const NT &lambda_prev) { return intersect_double_line_Vpoly(V, r, v, row, colno); } @@ -374,14 +371,14 @@ class VPolytope{ return std::pair (intersect_line_Vpoly(V, r, v, conv_comb, row, colno, false, false), 1); } - std::pair line_positive_intersect(const Point &r, const Point &v, const std::vector &Ar, - const std::vector &Av) { + std::pair line_positive_intersect(const Point &r, const Point &v, const VT &Ar, + const VT &Av) { return line_positive_intersect(r, v); } - std::pair line_positive_intersect(const Point &r, const Point &v, const std::vector &Ar, - const std::vector &Av, const NT &lambda_prev) { + std::pair line_positive_intersect(const Point &r, const Point &v, const VT &Ar, + const VT &Av, const NT &lambda_prev) { return line_positive_intersect(r, v);//, Ar, Av); } @@ -390,7 +387,7 @@ class VPolytope{ // with the V-polytope std::pair line_intersect_coord(const Point &r, const unsigned int rand_coord, - const std::vector &lamdas) { + const VT &lamdas) { std::vector temp(_d); temp[rand_coord]=1.0; @@ -405,7 +402,7 @@ class VPolytope{ const Point &r_prev, const unsigned int rand_coord, const unsigned int rand_coord_prev, - const std::vector &lamdas) { + const VT &lamdas) { return line_intersect_coord(r, rand_coord, lamdas); } @@ -441,14 +438,10 @@ class VPolytope{ return false; } unsigned int j; - std::vector temp(_d,NT(0)); + typename std::vector::iterator pointIt; for (int i=0; i 1.0) a = -a; - a = a/a.norm(); - - Point s(_d, std::vector(&a[0], a.data()+a.cols()*a.rows())); + a /= a.norm(); - s = ((-2.0 * v.dot(s)) * s); - v = s + v; + // compute reflection + a *= (-2.0 * v.dot(a)); + v += a; } void free_them_all() { diff --git a/include/convex_bodies/zonoIntersecthpoly.h b/include/convex_bodies/zonoIntersecthpoly.h index c39182efb..8ebae3148 100644 --- a/include/convex_bodies/zonoIntersecthpoly.h +++ b/include/convex_bodies/zonoIntersecthpoly.h @@ -3,6 +3,8 @@ // 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 ZONOINTERSECTHPOLY_H #define ZONOINTERSECTHPOLY_H @@ -63,16 +65,16 @@ class ZonoIntersectHPoly { std::max(polypair.second, zonopair.second)); } - std::pair line_intersect(Point &r, Point &v, std::vector &Ar, - std::vector &Av) { + std::pair line_intersect(Point &r, Point &v, VT &Ar, + VT &Av) { std::pair polypair = HP.line_intersect(r, v, Ar, Av); std::pair zonopair = Z.line_intersect(r, v); return std::pair(std::min(polypair.first, zonopair.first), std::max(polypair.second, zonopair.second)); } - std::pair line_intersect(Point &r, Point &v, std::vector &Ar, - std::vector &Av, NT &lambda_prev) { + std::pair line_intersect(Point &r, Point &v, VT &Ar, + VT &Av, NT &lambda_prev) { std::pair polypair = HP.line_intersect(r, v, Ar, Av, lambda_prev); std::pair zonopair = Z.line_intersect(r, v); return std::pair(std::min(polypair.first, zonopair.first), @@ -81,7 +83,7 @@ class ZonoIntersectHPoly { //First coordinate ray shooting intersecting convex body std::pair line_intersect_coord(Point &r,const unsigned int &rand_coord, - std::vector &lamdas) { + VT &lamdas) { std::pair polypair = HP.line_intersect_coord(r, rand_coord, lamdas); std::pair zonopair = Z.line_intersect_coord(r, rand_coord, lamdas); @@ -90,8 +92,8 @@ class ZonoIntersectHPoly { } - std::pair line_positive_intersect(Point &r, Point &v, std::vector &Ar, - std::vector &Av) { + std::pair line_positive_intersect(Point &r, Point &v, VT &Ar, + VT &Av) { std::pair polypair = HP.line_positive_intersect(r, v, Ar, Av); std::pair zonopair = Z.line_positive_intersect(r, v, Ar, Av); @@ -102,8 +104,8 @@ class ZonoIntersectHPoly { return std::pair(std::min(polypair.first, zonopair.first), facet); } - std::pair line_positive_intersect(Point &r, Point &v, std::vector &Ar, - std::vector &Av, NT &lambda_prev) { + std::pair line_positive_intersect(Point &r, Point &v, VT &Ar, + VT &Av, NT &lambda_prev) { std::pair polypair = HP.line_positive_intersect(r, v, Ar, Av, lambda_prev); std::pair zonopair = Z.line_positive_intersect(r, v, Ar, Av); int facet = HP.num_of_hyperplanes()+1; @@ -118,7 +120,7 @@ class ZonoIntersectHPoly { const Point &r_prev, const unsigned int &rand_coord, const unsigned int &rand_coord_prev, - std::vector &lamdas) { + VT &lamdas) { std::pair polypair = HP.line_intersect_coord(r, r_prev, rand_coord, rand_coord_prev, lamdas); std::pair zonopair = Z.line_intersect_coord(r, rand_coord, lamdas); diff --git a/include/convex_bodies/zpolytope.h b/include/convex_bodies/zpolytope.h index 7997eafdb..e79f5a122 100644 --- a/include/convex_bodies/zpolytope.h +++ b/include/convex_bodies/zpolytope.h @@ -4,6 +4,7 @@ // Copyright (c) 2018 Apostolos Chalkis //Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 program. +//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. // Licensed under GNU LGPL.3, see LICENCE file @@ -291,27 +292,27 @@ class Zonotope { // compute intersection point of ray starting from r and pointing to v // with the Zonotope - std::pair line_intersect(const Point &r, const Point &v, const std::vector &Ar, - const std::vector &Av) { + std::pair line_intersect(const Point &r, const Point &v, const VT& Ar, + const VT& Av) { return intersect_line_zono(V, r, v, conv_comb, colno); } // compute intersection point of ray starting from r and pointing to v // with the Zonotope - std::pair line_intersect(const Point &r, const Point &v, const std::vector &Ar, - const std::vector &Av, const NT &lambda_prev) { + std::pair line_intersect(const Point &r, const Point &v, const VT &Ar, + const VT &Av, const NT &lambda_prev) { return intersect_line_zono(V, r, v, conv_comb, colno); } - std::pair line_positive_intersect(const Point &r, const Point &v, const std::vector &Ar, - const std::vector &Av) { + std::pair line_positive_intersect(const Point &r, const Point &v, const VT &Ar, + const VT &Av) { return std::pair (intersect_line_Vpoly(V, r, v, conv_comb, row, colno, false, true), 1); } - std::pair line_positive_intersect(const Point &r, const Point &v, const std::vector &Ar, - const std::vector &Av, const NT &lambda_prev) { + std::pair line_positive_intersect(const Point &r, const Point &v, const VT &Ar, + const VT &Av, const NT &lambda_prev) { return line_positive_intersect(r, v, Ar, Av); } @@ -320,7 +321,7 @@ class Zonotope { // with the Zonotope std::pair line_intersect_coord(const Point &r, const unsigned int rand_coord, - const std::vector &lamdas) { + const VT &lamdas) { std::vector temp(_d,0); temp[rand_coord]=1.0; @@ -337,7 +338,7 @@ class Zonotope { const Point &r_prev, const unsigned int rand_coord, const unsigned int rand_coord_prev, - const std::vector &lamdas) { + const VT &lamdas) { return line_intersect_coord(r, rand_coord, lamdas); } @@ -386,17 +387,14 @@ class Zonotope { } VT a = Fmat.fullPivLu().kernel(); - NT sum = 0.0; - for (int k = 0; k < _d; ++k) sum += a(k)*p[k]; - if(sum<0.0) a = -1.0*a; + if(p.getCoefficients().dot(a) < 0.0) a = -1.0*a; a = a/a.norm(); - Point s(_d, std::vector(&a[0], a.data()+a.cols()*a.rows())); - - s = ((-2.0 * v.dot(s)) * s); - v = s + v; + // compute reflection + a *= (-2.0 * v.dot(a)); + v += a; } void free_them_all() { diff --git a/include/generators/h_polytopes_gen.h b/include/generators/h_polytopes_gen.h index f0e273d64..b0319860a 100644 --- a/include/generators/h_polytopes_gen.h +++ b/include/generators/h_polytopes_gen.h @@ -34,13 +34,8 @@ Polytope random_hpoly(unsigned int dim, unsigned int m, double seed = std::numer for(unsigned int i=0; i(dim); - pit = p.iter_begin(); - j = 0; - for ( ; pit!=p.iter_end(); ++pit, ++j) { - A(i,j) = *pit; - } + A.row(i) = p.getCoefficients(); b(i) = 10.0; - } Polytope HP; HP.init(dim, A, b); diff --git a/include/generators/known_polytope_generators.h b/include/generators/known_polytope_generators.h index df52215ba..f68060238 100644 --- a/include/generators/known_polytope_generators.h +++ b/include/generators/known_polytope_generators.h @@ -306,6 +306,7 @@ Polytope gen_skinny_cube(const unsigned int &dim, bool Vpoly = false) { return P; } + /* * ToDo: brkhoff polytope generator template diff --git a/include/lp_oracles/solve_lp.h b/include/lp_oracles/solve_lp.h index 094e777e2..a5a079afb 100644 --- a/include/lp_oracles/solve_lp.h +++ b/include/lp_oracles/solve_lp.h @@ -3,6 +3,7 @@ // Copyright (c) 2018 Vissarion Fisikopoulos, Apostolos Chalkis //Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 program. +//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. // VolEsti is free software: you can redistribute it and/or modify it // under the terms of the GNU Lesser General Public License as published by @@ -242,11 +243,15 @@ Point PointInIntersection(MT V1, MT V2, Point direction, bool &empty) { set_add_rowmode(lp, FALSE); /* rowmode should be turned off again when done building the model */ + const NT* direction_data = direction.getCoefficients().data(); + REAL* row_temp = row; + // set the bounds - typename std::vector::iterator pit = direction.iter_begin(); - for(j=0; j #include #include diff --git a/include/lp_oracles/zpolyoracles.h b/include/lp_oracles/zpolyoracles.h index 3b8e853ef..06ad2c686 100644 --- a/include/lp_oracles/zpolyoracles.h +++ b/include/lp_oracles/zpolyoracles.h @@ -1,6 +1,7 @@ #ifndef ZPOLYORACLES_H #define ZPOLYORACLES_H + #include #include #include diff --git a/include/rotating.h b/include/rotating.h index a29c705cb..366b8c5f1 100644 --- a/include/rotating.h +++ b/include/rotating.h @@ -3,6 +3,7 @@ // Copyright (c) 20012-2019 Vissarion Fisikopoulos // Copyright (c) 2019 Apostolos Chalkis + // Licensed under GNU LGPL.3, see LICENCE file #ifndef ROTATING_H diff --git a/include/rounding.h b/include/rounding.h index 90e5b8f5d..63cd369f6 100644 --- a/include/rounding.h +++ b/include/rounding.h @@ -4,6 +4,7 @@ // Copyright (c) 2018 Apostolos Chalkis //Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 program. +//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. // The functions in this header file call Bojan Nikolic implementation // of Todd and Yildirim algorithm in "On Khachiyan's Algorithm for the Computation of Minimum Volume Enclosing Ellipsoids", 2005 @@ -26,6 +27,7 @@ std::pair rounding_min_ellipsoid(Polytope &P , const std::pair rounding_min_ellipsoid(Polytope &P , const std::pair Ap(n,randPoints.size()); typename std::list::iterator rpit=randPoints.begin(); - typename std::vector::iterator qit; + for ( ; rpit!=randPoints.end(); rpit++, j++) { - qit = (*rpit).iter_begin(); i=0; - for ( ; qit!=(*rpit).iter_end(); qit++, i++){ - Ap(i,j)=double(*qit); + for (i=0 ; idimension(); i++){ + Ap(i,j)=double((*rpit)[i]); } } - boost::numeric::ublas::matrix Q(n,n); + boost::numeric::ublas::matrix Q(n,n); //TODO: remove dependence on ublas and copy to eigen boost::numeric::ublas::vector c2(n); size_t w=1000; KhachiyanAlgo(Ap,0.01,w,Q,c2); // call Khachiyan algorithm @@ -106,18 +107,16 @@ void get_vpoly_center(Polytope &P) { typedef typename Polytope::VT VT; typedef typename Polytope::PolytopePoint Point; - unsigned int n = P.dimension(), i, j = 0; + unsigned int n = P.dimension(); std::list randPoints; //ds for storing rand points P.get_points_for_rounding(randPoints); boost::numeric::ublas::matrix Ap(n,randPoints.size()); typename std::list::iterator rpit=randPoints.begin(); - typename std::vector::iterator qit; - for ( ; rpit!=randPoints.end(); rpit++, j++) { - qit = (*rpit).iter_begin(); i=0; - for ( ; qit!=(*rpit).iter_end(); qit++, i++){ - Ap(i,j)=double(*qit); + for (int j=0 ; rpit!=randPoints.end(); rpit++, j++) { + for (int i=0 ; idimension(); i++){ + Ap(i,j)=double((*rpit)[i]); } } boost::numeric::ublas::matrix Q(n,n); diff --git a/include/samplers/gaussian_samplers.h b/include/samplers/gaussian_samplers.h index 25da7455a..607e240fe 100644 --- a/include/samplers/gaussian_samplers.h +++ b/include/samplers/gaussian_samplers.h @@ -3,6 +3,7 @@ // Copyright (c) 2018 Vissarion Fisikopoulos, Apostolos Chalkis //Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 program. +//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. // VolEsti is free software: you can redistribute it and/or modify it // under the terms of the GNU Lesser General Public License as published by @@ -59,7 +60,7 @@ NT get_max_coord(const NT &l, const NT &u, const NT &a_i) { // Pick a point from the distribution exp(-a_i||x||^2) on the chord template -void rand_exp_range(Point &lower, Point &upper, const NT &a_i, Point &p, Parameters const& var) { +void rand_exp_range(Point const &lower, Point const & upper, const NT &a_i, Point &p, Parameters const& var) { typedef typename Parameters::RNGType RNGType; NT r, r_val, fn; const NT tol = 0.00000001; @@ -142,14 +143,14 @@ NT rand_exp_range_coord(const NT &l, const NT &u, const NT &a_i, Parameters cons // compute the first coordinate point -template +template void gaussian_first_coord_point(Polytope &P, Point &p, // a point to start Point &p_prev, // previous point unsigned int &coord_prev, // previous coordinate ray unsigned int walk_len, // number of steps for the random walk const NT &a_i, - std::vector &lamdas, + VT &lamdas, Parameters const& var) { typedef typename Parameters::RNGType RNGType; unsigned int n = var.n, rand_coord; @@ -175,14 +176,14 @@ void gaussian_first_coord_point(Polytope &P, // Compute the next point with target distribution the gaussian -template +template void gaussian_next_point(Polytope &P, Point &p, // a point to start Point &p_prev, // previous point unsigned int &coord_prev, // previous coordinate ray const unsigned int walk_len, // number of steps for the random walk const NT &a_i, - std::vector &lamdas, + VT &lamdas, Parameters const& var) { typedef typename Parameters::RNGType RNGType; unsigned int n = var.n, rand_coord; @@ -223,7 +224,10 @@ void rand_gaussian_point_generator(Polytope &P, RNGType &rng2 = var.rng; boost::random::uniform_int_distribution<> uidist(0, n - 1); - std::vector lamdas(P.num_of_hyperplanes(), NT(0)); + typedef Eigen::Matrix VT; + VT lamdas; + lamdas.setZero(P.num_of_hyperplanes()); + unsigned int rand_coord = uidist(rng2), coord_prev, rand_coord_prev; NT ball_rad = var.delta; Point p_prev = p; @@ -282,14 +286,14 @@ void gaussian_hit_and_run(Point &p, // hit-and-run with orthogonal directions and update -template +template void gaussian_hit_and_run_coord_update(Point &p, Point &p_prev, Polytope &P, unsigned int rand_coord, unsigned int rand_coord_prev, const NT &a_i, - std::vector &lamdas, + VT &lamdas, Parameters const& var) { std::pair bpair = P.line_intersect_coord(p, p_prev, rand_coord, rand_coord_prev, lamdas); NT dis = rand_exp_range_coord(p[rand_coord] + bpair.second, p[rand_coord] + bpair.first, a_i, var); @@ -310,7 +314,7 @@ void gaussian_ball_walk(Point & p, unsigned int n = P.dimension(); NT f_x, f_y, rnd; Point y = get_point_in_Dsphere(n, ball_rad); - y = y + p; + y += p; f_x = eval_exp(p, a_i); //unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); //RNGType rng(seed); diff --git a/include/samplers/samplers.h b/include/samplers/samplers.h index e0d555cf6..228b4c3cb 100644 --- a/include/samplers/samplers.h +++ b/include/samplers/samplers.h @@ -4,6 +4,7 @@ // Copyright (c) 2018 Apostolos Chalkis //Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 program. +//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. // Licensed under GNU LGPL.3, see LICENCE file @@ -11,26 +12,30 @@ #define RANDOM_SAMPLERS_H + + // Pick a random direction as a normilized vector template Point get_direction(const unsigned int dim) { boost::normal_distribution<> rdist(0,1); - std::vector Xs(dim,0); NT normal = NT(0); unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); RNGType rng(seed); + + Point p(dim); + NT* data = p.pointerToData(); + //RNGType rng2 = var.rng; - for (unsigned int i=0; i Point get_point_on_Dsphere(const unsigned int dim, const NT &radius){ Point p = get_direction(dim); - p = (radius == 0) ? p : radius * p; + if (radius != 0) p *= radius; return p; } @@ -55,7 +60,7 @@ Point get_point_in_Dsphere(const unsigned int dim, const NT &radius){ Point p = get_direction(dim); U = urdist(rng2); U = std::pow(U, 1.0/(NT(dim))); - p = (radius*U)*p; + p *= (radius*U); return p; } @@ -67,7 +72,7 @@ void ball_walk(Point &p, { //typedef typename Parameters::RNGType RNGType; Point y = get_point_in_Dsphere(p.dimension(), delta); - y = y + p; + y += p; if (P.is_in(y)==-1) p = y; } @@ -111,12 +116,18 @@ void boundary_rand_point_generator(Polytope &P, { typedef typename Parameters::RNGType RNGType; typedef typename Point::FT NT; + typedef typename Point::Coeff VT; + unsigned 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::vector lamdas(P.num_of_hyperplanes(), NT(0)), Av(P.num_of_hyperplanes(), NT(0)); + VT lamdas, Av; + + lamdas.setZero(P.num_of_hyperplanes()); + Av.setZero(P.num_of_hyperplanes()); + unsigned int rand_coord, rand_coord_prev; NT kapa, lambda; Point p_prev = p, p1(n), p2(n), v(n); @@ -132,7 +143,7 @@ void boundary_rand_point_generator(Polytope &P, v = get_direction(n); std::pair bpair = P.line_intersect(p, v, lamdas, Av); lambda = urdist(rng) * (bpair.first - bpair.second) + bpair.second; - p = (lambda * v) + p; + p += (lambda * v); } //hit_and_run(p, P, var); @@ -158,7 +169,7 @@ void boundary_rand_point_generator(Polytope &P, p1 = (bpair.first * v) + p; p2 = (bpair.second * v) + p; lambda = urdist(rng) * (bpair.first - bpair.second) + bpair.second; - p = (lambda * v) + p; + p += (lambda * v); } } @@ -184,7 +195,12 @@ void rand_point_generator(Polytope &P, boost::random::uniform_real_distribution<> urdist(0, 1); boost::random::uniform_int_distribution<> uidist(0, n - 1); - std::vector lamdas(P.num_of_hyperplanes(), NT(0)), Av(P.num_of_hyperplanes(), NT(0)); + typedef typename Point::Coeff VT; + VT lamdas, Av; + + lamdas.setZero(P.num_of_hyperplanes()); + Av.setZero(P.num_of_hyperplanes()); + unsigned int rand_coord, rand_coord_prev; NT kapa, ball_rad = var.delta, lambda; Point p_prev = p, v(n); @@ -201,7 +217,7 @@ void rand_point_generator(Polytope &P, v = get_direction(n); std::pair bpair = P.line_intersect(p, v, lamdas, Av); lambda = urdist(rng) * (bpair.first - bpair.second) + bpair.second; - p = (lambda * v) + p; + p += (lambda * v); //hit_and_run(p, P, var); } else { billiard_walk(P, p, var.diameter, lamdas, Av, lambda, var, true); @@ -221,7 +237,7 @@ void rand_point_generator(Polytope &P, v = get_direction(n); std::pair bpair = P.line_intersect(p, v, lamdas, Av, lambda); lambda = urdist(rng) * (bpair.first - bpair.second) + bpair.second; - p = (lambda * v) + p; + p += (lambda * v); //hit_and_run(p, P, var); } else { billiard_walk(P, p, var.diameter, lamdas, Av, lambda, var); @@ -250,7 +266,13 @@ void rand_point_generator(BallPoly &PBLarge, RNGType &rng = var.rng; boost::random::uniform_real_distribution<> urdist(0, 1); boost::random::uniform_int_distribution<> uidist(0, n - 1); - std::vector lamdas(PBLarge.num_of_hyperplanes(), NT(0)), Av(PBLarge.num_of_hyperplanes(), NT(0)); + + typedef typename Point::Coeff VT; + VT lamdas, Av; + + lamdas.setZero(PBLarge.num_of_hyperplanes()); + Av.setZero(PBLarge.num_of_hyperplanes()); + unsigned int rand_coord, rand_coord_prev; NT kapa, ball_rad = var.delta, lambda; Point p_prev = p, v(n); @@ -268,7 +290,7 @@ void rand_point_generator(BallPoly &PBLarge, v = get_direction(n); std::pair bpair = PBLarge.line_intersect(p, v, lamdas, Av); lambda = urdist(rng) * (bpair.first - bpair.second) + bpair.second; - p = (lambda * v) + p; + p += (lambda * v); //hit_and_run(p, PBLarge, var); } else { billiard_walk(PBLarge, p, var.diameter, lamdas, Av, lambda, var, true); @@ -287,7 +309,7 @@ void rand_point_generator(BallPoly &PBLarge, v = get_direction(n); std::pair bpair = PBLarge.line_intersect(p, v, lamdas, Av, lambda); lambda = urdist(rng) * (bpair.first - bpair.second) + bpair.second; - p = (lambda * v) + p; + p += (lambda * v); //hit_and_run(p, PBLarge, var); } else { billiard_walk(PBLarge, p, var.diameter, lamdas, Av, lambda, var); @@ -300,14 +322,14 @@ void rand_point_generator(BallPoly &PBLarge, } } -template +template void uniform_first_point(Polytope &P, Point &p, // a point to start Point &p_prev, // previous point unsigned int &coord_prev, // previous coordinate ray unsigned int walk_len, // number of steps for the random walk - std::vector &lamdas, - std::vector &Av, + VT &lamdas, + VT &Av, NT &lambda, const Parameters &var) { typedef typename Parameters::RNGType RNGType; @@ -332,7 +354,7 @@ void uniform_first_point(Polytope &P, v = get_direction(n); std::pair bpair = P.line_intersect(p, v, lamdas, Av); lambda = urdist(rng) * (bpair.first - bpair.second) + bpair.second; - p = (lambda * v) + p; + p += (lambda * v); } else { billiard_walk(P, p, var.diameter, lamdas, Av, lambda, var, true); } @@ -355,7 +377,7 @@ void uniform_first_point(Polytope &P, v = get_direction(n); std::pair bpair = P.line_intersect(p, v, lamdas, Av, lambda); lambda = urdist(rng) * (bpair.first - bpair.second) + bpair.second; - p = (lambda * v) + p; + p += (lambda * v); } } else { billiard_walk(P, p, var.diameter, lamdas, Av, lambda, var); @@ -364,14 +386,14 @@ void uniform_first_point(Polytope &P, -template +template void uniform_next_point(Polytope &P, Point &p, // a point to start Point &p_prev, // previous point unsigned int &coord_prev, // previous coordinate ray const unsigned int walk_len, // number of steps for the random walk - std::vector &lamdas, - std::vector &Av, + VT &lamdas, + VT &Av, NT &lambda, const Parameters &var) { typedef typename Parameters::RNGType RNGType; @@ -391,7 +413,7 @@ void uniform_next_point(Polytope &P, v = get_direction(n); std::pair bpair = P.line_intersect(p, v, lamdas, Av, lambda); lambda = urdist(rng) * (bpair.first - bpair.second) + bpair.second; - p = (lambda * v) + p; + p += (lambda * v); } } else if (var.bill_walk) { for (unsigned int j = 0; j < walk_len; j++) billiard_walk(P, p, var.diameter, lamdas, Av, lambda, var); @@ -427,22 +449,22 @@ void hit_and_run(Point &p, //hit-and-run with orthogonal directions and update -template +template void hit_and_run_coord_update(Point &p, Point &p_prev, Polytope &P, unsigned int rand_coord, unsigned int rand_coord_prev, const NT &kapa, - std::vector &lamdas) { + VT &lamdas) { std::pair bpair = P.line_intersect_coord(p, p_prev, rand_coord, rand_coord_prev, lamdas); p_prev = p; p.set_coord(rand_coord, p[rand_coord] + bpair.first + kapa * (bpair.second - bpair.first)); } -template -void billiard_walk(ConvexBody &P, Point &p, NT diameter, std::vector &Ar, std::vector &Av, NT &lambda_prev, +template +void billiard_walk(ConvexBody &P, Point &p, NT diameter, VT &Ar, VT &Av, NT &lambda_prev, Parameters &var, bool first = false) { typedef typename Parameters::RNGType RNGType; @@ -463,7 +485,7 @@ void billiard_walk(ConvexBody &P, Point &p, NT diameter, std::vector &Ar, st return; } lambda_prev = dl * pbpair.first; - p = (lambda_prev * v) + p; + p += (lambda_prev * v); T -= lambda_prev; P.compute_reflection(v, p, pbpair.second); } @@ -478,7 +500,7 @@ void billiard_walk(ConvexBody &P, Point &p, NT diameter, std::vector &Ar, st } lambda_prev = dl * pbpair.first; - p = (lambda_prev * v) + p; + p += (lambda_prev * v); T -= lambda_prev; P.compute_reflection(v, p, pbpair.second); it++; diff --git a/include/volume/cooling_balls.h b/include/volume/cooling_balls.h index 24436ceb8..399648070 100644 --- a/include/volume/cooling_balls.h +++ b/include/volume/cooling_balls.h @@ -3,6 +3,8 @@ // Copyright (c) 20012-2019 Vissarion Fisikopoulos // Copyright (c) 2018-2019 Apostolos Chalkis +//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. + #ifndef COOLING_BALLS_H #define COOLING_BALLS_H @@ -65,8 +67,9 @@ NT vol_cooling_balls(Polytope &P, UParameters &var, AParameters &var_ban, std::p // Save the radius of the Chebychev ball var.che_rad = radius; // Move the chebychev center to the origin and apply the same shifting to the polytope - VT c_e = Eigen::Map(&c.get_coeffs()[0], c.dimension()); - P.shift(c_e); + + P.shift(c.getCoefficients()); + if ( !get_sequence_of_polyballs(P, BallSet, ratios, N * nu, nu, lb, ub, radius, alpha, var, rmax) ){ return -1.0; diff --git a/include/volume/cooling_hpoly.h b/include/volume/cooling_hpoly.h index 57777bd5a..7b0025900 100644 --- a/include/volume/cooling_hpoly.h +++ b/include/volume/cooling_hpoly.h @@ -3,6 +3,7 @@ // Copyright (c) 20012-2018 Vissarion Fisikopoulos // Copyright (c) 2018 Apostolos Chalkis + #ifndef COOLING_HPOLY_H #define COOLING_HPOLY_H diff --git a/include/volume/exact_vols.h b/include/volume/exact_vols.h index 71525f63d..ceb97e18b 100644 --- a/include/volume/exact_vols.h +++ b/include/volume/exact_vols.h @@ -3,6 +3,8 @@ // 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. + // Licensed under GNU LGPL.3, see LICENCE file #ifndef ZONOTOPE_EXACT_VOL_H diff --git a/include/volume/volume.h b/include/volume/volume.h index f5aeac3a3..72c49aa79 100644 --- a/include/volume/volume.h +++ b/include/volume/volume.h @@ -5,6 +5,7 @@ // Copyright (c) 2018 Vissarion Fisikopoulos, Apostolos Chalkis //Contributed and/or modified by Apostolos Chalkis, as part of Google Summer of Code 2018 program. +//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. // Licensed under GNU LGPL.3, see LICENCE file @@ -78,8 +79,7 @@ NT volume(Polytope &P, } // Move the chebychev center to the origin and apply the same shifting to the polytope - VT c_e = Eigen::Map(&c.get_coeffs()[0], c.dimension()); - P.shift(c_e); + P.shift(c.getCoefficients()); c=Point(n); rnum=rnum/n_threads; @@ -288,8 +288,7 @@ NT volume_gaussian_annealing(Polytope &P, var.che_rad = radius; // Move the chebychev center to the origin and apply the same shifting to the polytope - VT c_e = Eigen::Map(&c.get_coeffs()[0], c.dimension()); - P.shift(c_e); + P.shift(c.getCoefficients()); // Initialization for the schedule annealing std::vector a_vals; @@ -321,7 +320,9 @@ NT volume_gaussian_annealing(Polytope &P, // Initialization for the approximation of the ratios unsigned int W = var.W, coord_prev, i=0; - std::vector last_W2(W,0), fn(mm,0), its(mm,0), lamdas(m,0); + std::vector last_W2(W,0), fn(mm,0), its(mm,0); + VT lamdas; + lamdas.setZero(m); vol=std::pow(M_PI/a_vals[0], (NT(n))/2.0)*std::abs(round_value); Point p(n), p_prev(n); // The origin is the Chebychev center of the Polytope viterator fnIt = fn.begin(), itsIt = its.begin(), avalsIt = a_vals.begin(), minmaxIt; diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index a5628e4f5..43fcef4ae 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -33,7 +33,7 @@ else () set(CMAKE_EXPORT_COMPILE_COMMANDS "ON") - #include_directories (BEFORE ../external/Eigen) + include_directories (BEFORE ../external/Eigen) include_directories (BEFORE ../external) include_directories (BEFORE ../external/minimum_ellipsoid) #include_directories (BEFORE ../include/cartesian_geom) @@ -50,7 +50,13 @@ else () include_directories (BEFORE ../include/lp_oracles) include_directories (BEFORE ../include/misc) - + #for Eigen + if (${CMAKE_VERSION} VERSION_LESS "3.12.0") + add_compile_options(-D "EIGEN_NO_DEBUG") + else () + add_compile_definitions("EIGEN_NO_DEBUG") + endif () + add_definitions(${CMAKE_CXX_FLAGS} "-std=c++11") # enable C++11 standard add_definitions(${CMAKE_CXX_FLAGS} "-O3") # optimization of the compiler #add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lgsl") @@ -68,6 +74,7 @@ else () add_library(test_main OBJECT test_main.cpp) add_executable (volume_test volume_test.cpp $) + add_executable (cheb_test chebychev_test.cpp $) #add_executable (rounding_test rounding_test.cpp $) add_executable (volumeCG_test volumeCG_test.cpp $)