Gaussian HMC -> length parameter integration in cooling gaussians - d…

…elta update on each step

examples/volume_cooling_gaussians_hmc/CMakeLists.txt

@@ -74,21 +74,8 @@ else ()
   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/lp_oracles)
-  include_directories (BEFORE ../../include/nlp_oracles)
-  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/misc)
-  include_directories (BEFORE ../../include/optimization)
-  include_directories (BEFORE ../../test)
+
 
   # for Eigen
   if (${CMAKE_VERSION} VERSION_LESS "3.12.0")
@@ -98,7 +85,6 @@ else ()
   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")

examples/volume_cooling_gaussians_hmc/volume_cooling_gaussians_hmc.cpp

@@ -1,25 +1,25 @@
 #include "Eigen/Eigen"
 #include <vector>
 #include "cartesian_geom/cartesian_kernel.h"
-#include "hpolytope.h"
-#include "known_polytope_generators.h"
+#include "convex_bodies/hpolytope.h"
+#include "generators/known_polytope_generators.h"
 #include "random_walks/random_walks.hpp"
-#include "volume_sequence_of_balls.hpp"
-#include "volume_cooling_gaussians.hpp"
-#include "volume_cooling_balls.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.h"
+#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 calculateAndVerifyVolume(HPOLYTOPE& polytope, const std::string& description, NT expectedVolume, NT epsilon, NT L) {
-    int walk_len = 200;
+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";
@@ -40,19 +40,19 @@ int main() {
 
     // 3-dimensional cube
     HPOLYTOPE cube = generate_cube<HPOLYTOPE>(3, false);
-    calculateAndVerifyVolume(cube, "3-dimensional cube", 1.0, epsilon, 10);
+    verify_volume(cube, "3-dimensional cube", 1.0, epsilon, 4);
 
     // 3-dimensional cross polytope
     HPOLYTOPE crossPolytope = generate_cross<HPOLYTOPE>(3, false);
-    calculateAndVerifyVolume(crossPolytope, "3-dimensional cross polytope", 4.0 / 6.0, epsilon, 10);
+    verify_volume(crossPolytope, "3-dimensional cross polytope", 4.0 / 6.0, epsilon, 4);
 
     // 3-dimensional simplex
     HPOLYTOPE simplex = generate_simplex<HPOLYTOPE>(3, false);
-    calculateAndVerifyVolume(simplex, "3-dimensional simplex", 1.0 / 6.0, epsilon, 4);
+    verify_volume(simplex, "3-dimensional simplex", 1.0 / 6.0, epsilon, 10);
 
-    // birkhoff polytope dimension 3 
+    // 4-dimensional birkhoff  
     HPOLYTOPE birkhoffPolytope = generate_birkhoff<HPOLYTOPE>(3);
-    calculateAndVerifyVolume(birkhoffPolytope, "Birkhoff polytope (dim 3)", -1, epsilon, 8);  // Theoretical volume not easily available
+    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,8 +224,7 @@ private :
         }
     }
 
-    inline void update_position(Point& p, Point& v, const NT& T, const NT& omega) {
-
+    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);

include/volume/volume_cooling_gaussians.hpp

@@ -45,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
 
@@ -178,7 +188,7 @@ NT get_next_gaussian(Polytope& P,
 
 template
 <
-    typename WalkTypePolicy,
+    typename WalkType,
     typename RandomPointGenerator,
     typename Polytope,
     typename NT,
@@ -194,7 +204,7 @@ void compute_annealing_schedule(Polytope& P,
                                 NT const& error,
                                 std::vector<NT>& a_vals,
                                 RandomNumberGenerator& rng,
-                                typename std::conditional<std::is_same<WalkTypePolicy, GaussianHamiltonianMonteCarloExactWalk>::value,
+                                typename std::conditional<std::is_same<WalkType, GaussianHamiltonianMonteCarloExactWalk>::value,
                                 typename GaussianHamiltonianMonteCarloExactWalk::parameters, void*>::type walk_params = nullptr)
 {
     typedef typename Polytope::PointType Point;
@@ -228,16 +238,27 @@ void compute_annealing_schedule(Polytope& P,
         NT curr_fn = 0;
         NT curr_its = 0;
 
-        if constexpr (std::is_same<WalkTypePolicy, GaussianHamiltonianMonteCarloExactWalk>::value) {
+        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)));
+
+
 
             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 {
-            WalkTypePolicy walk(P, p, a_vals[it], rng);
+
+            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)));
+
 
             for (unsigned int j = 0; j < totalSteps; ++j) {
                 walk.template apply(P, p, a_vals[it], walk_length, rng);
@@ -261,8 +282,6 @@ void compute_annealing_schedule(Polytope& P,
     }
 }
 
-
-
 template <typename NT>
 struct gaussian_annealing_parameters
 {
@@ -336,7 +355,6 @@ double volume_cooling_gaussians(Polytope& Pin,
     unsigned int N = parameters.N;
 
     // Construct the necesary parameters if this is the case
-    // Construct the necessary 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);
@@ -346,7 +364,8 @@ double volume_cooling_gaussians(Polytope& Pin,
             compute_annealing_schedule<WalkType, RandomPointGenerator>(
                 P, ratio, C, parameters.frac, N, walk_length, radius, error, a_vals, rng);
         }
-        } else {
+    } 
+    else {
             compute_annealing_schedule<WalkType, RandomPointGenerator>(
                 P, ratio, C, parameters.frac, N, walk_length, radius, error, a_vals, rng);
     }