Gaussian HMC - integration within cooling gaussians and complexity improvements

examples/volume_cooling_gaussians_hmc/CMakeLists.txt

@@ -0,0 +1,99 @@
+# VolEsti (volume computation and sampling library)
+# Copyright (c) 2012-2021 Vissarion Fisikopoulos
+# Copyright (c) 2018-2021 Apostolos Chalkis
+# Contributed and/or modified by Marios Papachristou, as part of Google Summer of Code 2020 program.
+# 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)
+
+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 ../../include)
+
+
+  # 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} "-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_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lgslcblas")
+  #add_definitions( "-O3 -lgsl -lm -ldl -lgslcblas" )
+
+  add_executable (volume_cooling_gaussians_hmc volume_cooling_gaussians_hmc.cpp)
+  target_link_libraries(volume_cooling_gaussians_hmc ${LP_SOLVE} ${IFOPT} ${IFOPT_IPOPT} ${GMP} ${MPSOLVE} ${PTHREAD} ${FFTW3})
+
+endif()

examples/volume_cooling_gaussians_hmc/volume_cooling_gaussians_hmc.cpp

@@ -0,0 +1,59 @@
+#include "Eigen/Eigen"
+#include <vector>
+#include "cartesian_geom/cartesian_kernel.h"
+#include "convex_bodies/hpolytope.h"
+#include "generators/known_polytope_generators.h"
+#include "random_walks/random_walks.hpp"
+#include "volume/volume_sequence_of_balls.hpp"
+#include "volume/volume_cooling_gaussians.hpp"
+#include "volume/volume_cooling_balls.hpp"
+
+#include <iostream>
+#include <fstream>
+#include "misc/misc.h"
+
+typedef double NT;
+typedef Cartesian<NT> Kernel;
+typedef typename Kernel::Point Point;
+typedef BoostRandomNumberGenerator<boost::mt19937, NT, 3> RNGType;
+typedef HPolytope<Point> HPOLYTOPE;
+
+void verify_volume(HPOLYTOPE& polytope, const std::string& description, NT expectedVolume, NT epsilon, NT L) {
+    int walk_len = 1;
+    NT e = 0.1;
+
+    std::cout << "Calculating volume for " << description << ":\n";
+    NT volume = volume_cooling_gaussians<GaussianHamiltonianMonteCarloExactWalk, RNGType>(polytope, e, walk_len, L);
+
+    std::cout << "Estimated volume: " << volume << "\n";
+    std::cout << "Expected volume: " << expectedVolume << "\n";
+    NT error = std::abs(volume - expectedVolume);
+
+    if (error <= epsilon)
+        std::cout << "Result is within the acceptable error margin (" << epsilon << ").\n\n";
+    else
+        std::cout << "Error (" << error << ") is larger than acceptable margin (" << epsilon << ").\n\n";
+}
+
+int main() {
+    const NT epsilon = 0.1;
+
+    // 3-dimensional cube
+    HPOLYTOPE cube = generate_cube<HPOLYTOPE>(3, false);
+    verify_volume(cube, "3-dimensional cube", 1.0, epsilon, 4);
+
+    // 3-dimensional cross polytope
+    HPOLYTOPE crossPolytope = generate_cross<HPOLYTOPE>(3, false);
+    verify_volume(crossPolytope, "3-dimensional cross polytope", 4.0 / 6.0, epsilon, 4);
+
+    // 3-dimensional simplex
+    HPOLYTOPE simplex = generate_simplex<HPOLYTOPE>(3, false);
+    verify_volume(simplex, "3-dimensional simplex", 1.0 / 6.0, epsilon, 10);
+
+    // 4-dimensional birkhoff  
+    HPOLYTOPE birkhoffPolytope = generate_birkhoff<HPOLYTOPE>(3);
+    verify_volume(birkhoffPolytope, "Birkhoff polytope (dim 4)", 1.0 / 16.0, epsilon, 2);  // Theoretical volume not easily available
+
+    return 0;
+}
+

include/random_walks/gaussian_hamiltonian_monte_carlo_exact_walk.hpp

@@ -224,22 +224,22 @@ private :
         }
     }
 
-    inline void update_position(Point &p, Point &v, NT const& T, NT const& omega)
-    {
-        NT C, Phi;
-        for (size_t i = 0; i < p.dimension(); i++)
-        {
-            C = std::sqrt(p[i] * p[i] + (v[i] * v[i]) / (omega * omega));
-            Phi = std::atan((-v[i]) / (p[i] * omega));
-            if (v[i] < 0.0 && Phi < 0.0) {
-                Phi += M_PI;
-            } else if (v[i] > 0.0 && Phi > 0.0) {
-                Phi -= M_PI;
-            }
-            p.set_coord(i, C * std::cos(omega * T + Phi));
-            v.set_coord(i, -C * omega * std::sin(omega * T + Phi));
-        }
+    inline void update_position(Point &p, Point &v, NT const& T, NT const& omega) {
+        NT next_p, next_v;
+
+        NT sinVal = std::sin(omega * T);
+        NT cosVal = std::cos(omega * T);
 
+        NT factor1 = sinVal / omega;
+        NT factor2 = -omega * sinVal;
+
+        for (size_t i = 0; i < p.dimension(); i++) {
+            next_p = cosVal * p[i] + v[i] * factor1;
+            next_v = factor2 * p[i] + cosVal * v[i];
+
+            p.set_coord(i, next_p);
+            v.set_coord(i, next_v);
+        }
     }
 
     inline double get_max_distance(std::vector<Point> &pointset, Point const& q, double &rad)
@@ -271,5 +271,4 @@ private :
 };
 
 
-#endif // RANDOM_WALKS_GAUSSIAN_HMC_WALK_HPP
-
+#endif // RANDOM_WALKS_GAUSSIAN_HMC_WALK_HPP

include/volume/volume_cooling_gaussians.hpp

@@ -20,8 +20,7 @@
 
 #include "cartesian_geom/cartesian_kernel.h"
 #include "random_walks/gaussian_helpers.hpp"
-#include "random_walks/gaussian_ball_walk.hpp"
-#include "random_walks/gaussian_cdhr_walk.hpp"
+#include "random_walks/random_walks.hpp"
 #include "sampling/random_point_generators.hpp"
 #include "volume/math_helpers.hpp"
 
@@ -46,6 +45,16 @@ struct update_delta<GaussianBallWalk::Walk<Polytope, RandomNumberGenerator>>
     }
 };
 
+template <typename Polytope, typename RandomNumberGenerator>
+struct update_delta<GaussianHamiltonianMonteCarloExactWalk::Walk<Polytope, RandomNumberGenerator>>
+{
+    template <typename NT>
+    static void apply(GaussianHamiltonianMonteCarloExactWalk::Walk<Polytope, RandomNumberGenerator>& walk, NT delta)
+    {
+        walk.update_delta(delta);
+    }
+};
+
 
 ////////////////////////////// Algorithms
 
@@ -177,7 +186,6 @@ NT get_next_gaussian(Polytope& P,
     return last_a * std::pow(ratio, k);
 }
 
-// Compute the sequence of spherical gaussians
 template
 <
     typename WalkType,
@@ -195,7 +203,9 @@ void compute_annealing_schedule(Polytope& P,
                                 NT const& chebychev_radius,
                                 NT const& error,
                                 std::vector<NT>& a_vals,
-                                RandomNumberGenerator& rng)
+                                RandomNumberGenerator& rng,
+                                typename std::conditional<std::is_same<WalkType, GaussianHamiltonianMonteCarloExactWalk>::value,
+                                typename GaussianHamiltonianMonteCarloExactWalk::parameters, void*>::type walk_params = nullptr)
 {
     typedef typename Polytope::PointType Point;
     typedef typename Polytope::VT VT;
@@ -204,56 +214,64 @@ void compute_annealing_schedule(Polytope& P,
     get_first_gaussian(P, frac, chebychev_radius, error, a_vals);
 
 #ifdef VOLESTI_DEBUG
-    std::cout<<"first gaussian computed\n"<<std::endl;
+    std::cout << "first gaussian computed\n" << std::endl;
 #endif
 
     NT a_stop = 0.0;
     const NT tol = 0.001;
     unsigned int it = 0;
     unsigned int n = P.dimension();
-    const unsigned int totalSteps = ((int)150/((1.0 - frac) * error))+1;
-
-    if (a_vals[0]<a_stop) a_vals[0] = a_stop;
-
-#ifdef VOLESTI_DEBUG
-    std::cout<<"Computing the sequence of gaussians..\n"<<std::endl;
-#endif
+    const unsigned int totalSteps = static_cast<unsigned int>(150 / ((1.0 - frac) * error)) + 1;
 
     Point p(n);
 
     while (true)
     {
         // Compute the next gaussian
         NT next_a = get_next_gaussian<RandomPointGenerator>
-                      (P, p, a_vals[it], N, ratio, C, walk_length, rng);
+                    (P, p, a_vals[it], N, ratio, C, walk_length, rng);
+
+#ifdef VOLESTI_DEBUG
+        std::cout << "Processing a_vals[" << it << "] = " << a_vals[it] << ", next_a = " << next_a << std::endl;
+#endif
 
         NT curr_fn = 0;
         NT curr_its = 0;
-        auto steps = totalSteps;
 
-        WalkType walk(P, p, a_vals[it], rng);
-        //TODO: test update delta here?
+        if constexpr (std::is_same<WalkType, GaussianHamiltonianMonteCarloExactWalk>::value) {
+            GaussianHamiltonianMonteCarloExactWalk::Walk<Polytope, RandomNumberGenerator> walk(P, p, a_vals[it], rng, walk_params);
+
+            update_delta<decltype(walk)>::apply(walk, 4.0 * chebychev_radius / std::sqrt(std::max(NT(1.0), a_vals[it]) * NT(n)));
 
-        update_delta<WalkType>
+
+
+            for (unsigned int j = 0; j < totalSteps; ++j) {
+                walk.template apply(P, p, a_vals[it], walk_length, rng);
+                curr_its++;
+                curr_fn += eval_exp(p, next_a) / eval_exp(p, a_vals[it]);
+            }
+        } else {
+
+            WalkType walk(P, p, a_vals[it], rng);
+            //TODO: test update delta here?
+
+            update_delta<WalkType>
                 ::apply(walk, 4.0 * chebychev_radius
                         / std::sqrt(std::max(NT(1.0), a_vals[it]) * NT(n)));
 
-        // Compute some ratios to decide if this is the last gaussian
-        for (unsigned  int j = 0; j < steps; j++)
-        {
-            walk.template apply(P, p, a_vals[it], walk_length, rng);
-            curr_its += 1.0;
-            curr_fn += eval_exp(p, next_a) / eval_exp(p, a_vals[it]);
+
+            for (unsigned int j = 0; j < totalSteps; ++j) {
+                walk.template apply(P, p, a_vals[it], walk_length, rng);
+                curr_its++;
+                curr_fn += eval_exp(p, next_a) / eval_exp(p, a_vals[it]);
+            }
         }
 
-        // Remove the last gaussian.
-        // Set the last a_i equal to zero
-        if (next_a>0 && curr_fn/curr_its>(1.0+tol))
-        {
+        // pick the last gaussian
+        if (next_a > 0 && curr_fn / curr_its > (1.0 + tol)) {
             a_vals.push_back(next_a);
             it++;
-        } else if (next_a <= 0)
-        {
+        } else if (next_a <= 0) {
             a_vals.push_back(a_stop);
             it++;
             break;
@@ -292,7 +310,8 @@ template
 double volume_cooling_gaussians(Polytope& Pin,
                                 RandomNumberGenerator& rng,
                                 double const& error = 0.1,
-                                unsigned int const& walk_length = 1)
+                                unsigned int const& walk_length = 1,
+                                double L = -1)
 {
     typedef typename Polytope::PointType Point;
     typedef typename Point::FT NT;
@@ -335,11 +354,22 @@ double volume_cooling_gaussians(Polytope& Pin,
     NT C = parameters.C;
     unsigned int N = parameters.N;
 
-    compute_annealing_schedule
-    <
-        WalkType,
-        RandomPointGenerator
-    >(P, ratio, C, parameters.frac, N, walk_length, radius, error, a_vals, rng);
+    // Construct the necesary parameters if this is the case
+    if constexpr (std::is_same<WalkTypePolicy, GaussianHamiltonianMonteCarloExactWalk>::value) {
+        if (L > 0) {
+            typename GaussianHamiltonianMonteCarloExactWalk::parameters walk_params(L, true, 0, false);
+            compute_annealing_schedule<GaussianHamiltonianMonteCarloExactWalk, RandomPointGenerator>(
+                P, ratio, C, parameters.frac, N, walk_length, radius, error, a_vals, rng, walk_params);
+        } else {
+            compute_annealing_schedule<WalkType, RandomPointGenerator>(
+                P, ratio, C, parameters.frac, N, walk_length, radius, error, a_vals, rng);
+        }
+    } 
+    else {
+            compute_annealing_schedule<WalkType, RandomPointGenerator>(
+                P, ratio, C, parameters.frac, N, walk_length, radius, error, a_vals, rng);
+    }
+
 
 #ifdef VOLESTI_DEBUG
     std::cout<<"All the variances of schedule_annealing computed in = "
@@ -487,4 +517,4 @@ double volume_cooling_gaussians(Polytope &Pin,
     return volume_cooling_gaussians<WalkTypePolicy>(Pin, rng, error, walk_length);
 }
 
-#endif // VOLUME_COOLING_GAUSSIANS_HPP
+#endif // VOLUME_COOLING_GAUSSIANS_HPP