Gaussian accelerated billiard walk for volume cooling ellipsoids algorithm.

examples/sampling-hpolytope-with-billiard-walks/.gitignore

@@ -0,0 +1,2 @@
+sampler
+*walk.txt

examples/sampling-hpolytope-with-billiard-walks/CMakeLists.txt

@@ -0,0 +1,115 @@
+# Licensed under GNU LGPL.3, see LICENCE file
+# VolEsti (volume computation and sampling library)
+
+# Copyright (c) 2012-2021 Vissarion Fisikopoulos
+# Copyright (c) 2018-2021 Apostolos Chalkis
+# Contributed and/or modified by Vaibhav Thakkar
+
+# Licensed under GNU LGPL.3, see LICENCE file
+
+project( VolEsti )
+
+
+CMAKE_MINIMUM_REQUIRED(VERSION 3.11)
+
+set(CMAKE_ALLOW_LOOSE_LOOP_CONSTRUCTS true)
+
+if(COMMAND cmake_policy)
+  cmake_policy(SET CMP0003 NEW)
+endif(COMMAND cmake_policy)
+
+
+option(DISABLE_NLP_ORACLES "Disable non-linear oracles (used in collocation)" ON)
+option(BUILTIN_EIGEN "Use eigen from ../../external" OFF)
+
+
+if(DISABLE_NLP_ORACLES)
+  add_definitions(-DDISABLE_NLP_ORACLES)
+else()
+  find_library(IFOPT NAMES libifopt_core.so PATHS /usr/local/lib)
+  find_library(IFOPT_IPOPT NAMES libifopt_ipopt.so PATHS /usr/local/lib)
+  find_library(GMP NAMES libgmp.so PATHS /usr/lib/x86_64-linux-gnu /usr/lib/i386-linux-gnu)
+  find_library(MPSOLVE NAMES libmps.so PATHS /usr/local/lib)
+  find_library(PTHREAD NAMES libpthread.so PATHS /usr/lib/x86_64-linux-gnu /usr/lib/i386-linux-gnu)
+  find_library(FFTW3 NAMES libfftw3.so.3 PATHS /usr/lib/x86_64-linux-gnu /usr/lib/i386-linux-gnu)
+
+  if (NOT IFOPT)
+
+    message(FATAL_ERROR "This program requires the ifopt library, and will not be compiled.")
+
+  elseif (NOT GMP)
+
+    message(FATAL_ERROR "This program requires the gmp library, and will not be compiled.")
+
+  elseif (NOT MPSOLVE)
+
+    message(FATAL_ERROR "This program requires the mpsolve library, and will not be compiled.")
+
+  elseif (NOT FFTW3)
+
+    message(FATAL_ERROR "This program requires the fftw3 library, and will not be compiled.")
+
+  else()
+    message(STATUS "Library ifopt found: ${IFOPT}")
+    message(STATUS "Library gmp found: ${GMP}")
+    message(STATUS "Library mpsolve found: ${MPSOLVE}")
+    message(STATUS "Library fftw3 found:" ${FFTW3})
+
+  endif(NOT IFOPT)
+
+endif(DISABLE_NLP_ORACLES)
+
+include("../../external/cmake-files/Eigen.cmake")
+GetEigen()
+
+include("../../external/cmake-files/Boost.cmake")
+GetBoost()
+
+include("../../external/cmake-files/LPSolve.cmake")
+GetLPSolve()
+
+# Find lpsolve library
+find_library(LP_SOLVE NAMES liblpsolve55.so PATHS /usr/lib/lp_solve)
+
+if (NOT LP_SOLVE)
+  message(FATAL_ERROR "This program requires the lp_solve library, and will not be compiled.")
+else ()
+  message(STATUS "Library lp_solve found: ${LP_SOLVE}")
+
+  set(CMAKE_EXPORT_COMPILE_COMMANDS "ON")
+
+  include_directories (BEFORE ../../external)
+  include_directories (BEFORE ../../external/minimum_ellipsoid)
+  include_directories (BEFORE ../../include/generators)
+  include_directories (BEFORE ../../include/volume)
+  include_directories (BEFORE ../../include)
+  include_directories (BEFORE ../../include/convex_bodies)
+  include_directories (BEFORE ../../include/random_walks)
+  include_directories (BEFORE ../../include/annealing)
+  include_directories (BEFORE ../../include/ode_solvers)
+  include_directories (BEFORE ../../include/root_finders)
+  include_directories (BEFORE ../../include/samplers)
+  include_directories (BEFORE ../../include/lp_oracles)
+  include_directories (BEFORE ../../include/nlp_oracles)
+  include_directories (BEFORE ../../include/misc)
+  include_directories (BEFORE ../../include/optimization)
+
+  # 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")
+  add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lm")
+  add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-ldl")
+  add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-DBOOST_NO_AUTO_PTR")
+
+  add_executable (sampler sampler.cpp)
+  TARGET_LINK_LIBRARIES(sampler ${LP_SOLVE})
+
+endif()

examples/sampling-hpolytope-with-billiard-walks/README.md

@@ -0,0 +1,30 @@
+## Compilation
+Build the example by running the following commands in this directory.
+
+```bash
+cmake . -DLP_SOLVE=_PATH_TO_LIB_FILE
+make
+```  
+You have to specify the path to liblpsolve55.so/dll/dylib.  
+For example: -DLP_SOLVE=/usr/lib/lpsolve/liblpsolve55.so
+
+## Usage:
+```bash
+ ./sampler
+```
+After this, you can run `python3 plot.py` to plot the sampled points.
+
+## Sample instance:
+Currently there is only sample instance (number of sample points = 1000) in the form `Ax<=b`, representing the following inequalities
+```
+x >= -1
+y >= -1
+x-y <= 8
+x-y >= -8
+x+y <= 50
+```
+
+**Sampled points:**
+![sampled_points_1](images/uniform_bw.png)
+![sampled_points_2](images/uniform_acc_bw.png)
+![sampled_points_3](images/gaussian_bw.png)

examples/sampling-hpolytope-with-billiard-walks/plot.py

@@ -0,0 +1,43 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+A = np.array([[-1, 0],
+              [0, -1],
+              [1, -1],
+              [-1, 1],
+              [ 1, 1],
+            ])
+b = np.array([1,
+              1,
+              8,
+              8,
+              50])
+
+
+def get_polytope_vals(x1, x2):
+    X = np.vstack([x1, x2]).T
+    vals = [np.max((A @ x) - b) for x in X]
+    return np.array(vals)
+
+
+# plot data
+filenames = ["uniform_billiard_walk.txt", "uniform_accelerated_billiard_walk.txt", "gaussian_billiard_walk.txt"]
+for filename in filenames:
+    fig = plt.figure()
+
+    # plot polytope
+    x1_vec, x2_vec = np.meshgrid(np.linspace(-1, 30, 100), np.linspace(-1, 30, 100))
+    vals = get_polytope_vals(x1_vec.ravel(), x2_vec.ravel())
+    vals = vals.reshape(x1_vec.shape)
+    plt.contourf(x1_vec, x2_vec, vals, levels=[-100, 0], cmap="Greys_r", alpha=0.3)
+
+    # plot points
+    data = np.genfromtxt(filename, delimiter=' ')
+    plt.plot(data[:, 0], data[:, 1], 'x')
+    plt.xlim(-2, 30)
+    plt.ylim(-2, 30)
+    plt.title(filename[:-4])
+
+# show all plots
+plt.show() 

examples/sampling-hpolytope-with-billiard-walks/sampler.cpp

@@ -0,0 +1,121 @@
+#include <iostream>
+#include <vector>
+#include <fstream>
+#include <boost/random.hpp>
+#include "Eigen/Eigen"
+#include "cartesian_geom/cartesian_kernel.h"
+#include "sampling/random_point_generators.hpp"
+#include "random_walks/random_walks.hpp"
+#include "preprocess/max_inscribed_ellipsoid.hpp"
+#include "convex_bodies/ellipsoid.h"
+#include "convex_bodies/hpolytope.h"
+
+
+typedef double NT;
+typedef Cartesian <NT> Kernel;
+typedef typename Kernel::Point Point;
+typedef Eigen::Matrix<NT, Eigen::Dynamic, Eigen::Dynamic> MT;
+typedef Eigen::Matrix<NT, Eigen::Dynamic, 1> VT;
+typedef BoostRandomNumberGenerator<boost::mt19937, NT> RNGType;
+typedef HPolytope <Point> HPOLYTOPE;
+
+// walk policy
+PushBackWalkPolicy push_back_policy;
+
+// different billiard walks
+typedef BilliardWalk::template Walk<HPOLYTOPE, RNGType> BilliardWalkType;
+typedef AcceleratedBilliardWalk::template Walk<HPOLYTOPE, RNGType> AcceleratedBilliardWalkType;
+typedef GaussianAcceleratedBilliardWalk::template Walk<HPOLYTOPE, RNGType> GaussianAcceleratedWalkType;
+
+
+void write_to_file(std::string filename, std::vector<Point> const& randPoints) {
+    std::ofstream out(filename);
+    auto coutbuf = std::cout.rdbuf(out.rdbuf()); //save and redirect
+    for(int i=0; i<randPoints.size(); ++i)
+        randPoints[i].print();
+    std::cout.rdbuf(coutbuf); //reset to standard output again
+}
+
+
+void sample_using_uniform_billiard_walk(HPOLYTOPE const& HP, RNGType& rng, unsigned int walk_len, unsigned int num_points) {
+    std::string filename = "uniform_billiard_walk.txt";
+    typedef RandomPointGenerator<BilliardWalkType> Generator;
+
+    std::vector<Point> randPoints;
+    Point q(HP.dimension()); // origin
+    Generator::apply(HP, q, num_points, walk_len,
+                     randPoints, push_back_policy, rng);
+    write_to_file(filename, randPoints);
+}
+
+
+void sample_using_uniform_accelerated_billiard_walk(HPOLYTOPE const& HP, RNGType& rng, unsigned int walk_len, unsigned int num_points) {
+    std::string filename = "uniform_accelerated_billiard_walk.txt";
+    typedef RandomPointGenerator<AcceleratedBilliardWalkType> Generator;
+
+    std::vector<Point> randPoints;
+    Point q(HP.dimension()); // origin
+    Generator::apply(HP, q, num_points, walk_len,
+                     randPoints, push_back_policy, rng);
+    write_to_file(filename, randPoints);
+}
+
+
+void sample_using_gaussian_billiard_walk(HPOLYTOPE const& HP, RNGType& rng, unsigned int walk_len, unsigned int num_points) {
+    std::string filename = "gaussian_billiard_walk.txt";
+    typedef MultivariateGaussianRandomPointGenerator<GaussianAcceleratedWalkType> Generator;
+
+    std::vector<Point> randPoints;
+    Point q(HP.dimension()); // origin
+
+    // ----------- Get inscribed ellipsoid --------------------------------
+    typedef Ellipsoid<Point> EllipsoidType;
+    unsigned int max_iter = 150;
+    NT tol = std::pow(10, -6.0), reg = std::pow(10, -4.0);
+    VT x0 = q.getCoefficients();
+    std::pair<std::pair<MT, VT>, bool> inscribed_ellipsoid_res = max_inscribed_ellipsoid<MT>(HP.get_mat(),
+                                                                                             HP.get_vec(),
+                                                                                             x0,
+                                                                                             max_iter,
+                                                                                             tol,
+                                                                                             reg);
+    if (!inscribed_ellipsoid_res.second) // not converged
+        throw std::runtime_error("max_inscribed_ellipsoid not converged");
+
+    MT A_ell = inscribed_ellipsoid_res.first.first.inverse();
+    EllipsoidType inscribed_ellipsoid(A_ell);
+    // --------------------------------------------------------------------
+
+    Generator::apply(HP, q, inscribed_ellipsoid, num_points, walk_len,
+                     randPoints, push_back_policy, rng);
+    write_to_file(filename, randPoints);
+}
+
+
+HPOLYTOPE visualization_example() {
+    MT Ab(5, 3);
+    Ab <<  -1,  0, 1,     // x >= -1
+            0, -1, 1,     // y >= -1
+            1, -1, 8,     // x-y <= 8
+           -1,  1, 8,     // x-y >= -8
+            1,  1, 50;    // x+y <= 50
+
+    MT A = Ab.block(0, 0, 5, 2);
+    VT b = Ab.col(2);
+    return HPOLYTOPE(2, A, b);
+}
+
+
+int main(int argc, char const *argv[]) {
+    unsigned int walk_len = 10;
+    unsigned int num_points = 1000;
+
+    HPOLYTOPE HP = visualization_example();
+    HP.ComputeInnerBall();
+
+    RNGType rng(HP.dimension());
+    sample_using_uniform_billiard_walk(HP, rng, walk_len, num_points);
+    sample_using_uniform_accelerated_billiard_walk(HP, rng, walk_len, num_points);
+    sample_using_gaussian_billiard_walk(HP, rng, walk_len, num_points);
+    return 0;
+}