CRHMC with Cooling Non-SphericalGaussians (#330)

* Initial CRHMC with Cooling Non-SphericalGaussians

* Improved CRHMC burnIn logic, with additional random hpoly examples

* Removed hessian related lines

examples/crhmc_cooling_nonspherical_gaussians/CMakeLists.txt

@@ -0,0 +1,132 @@
+# VolEsti (volume computation and sampling library)
+
+# Copyright (c) 2012-2024 Vissarion Fisikopoulos
+# Copyright (c) 2018-2024 Apostolos Chalkis
+# Copyright (c) 2024 Vladimir Necula
+
+# Contributed and/or modified by Vladimir Necula, as part of Google Summer of
+# Code 2024 program.
+
+# Licensed under GNU LGPL.3, see LICENCE file
+
+project( VolEsti )
+
+
+CMAKE_MINIMUM_REQUIRED(VERSION 3.11)
+
+set(CMAKE_ALLOW_LOOSE_LOOP_CONSTRUCTS true)
+
+# Locate Intel MKL root (in case it is enabled)
+if (APPLE)
+	set(MKLROOT /opt/intel/oneapi/mkl/latest)
+elseif(UNIX)
+  #set(MKLROOT /opt/intel/oneapi/mkl/latest)
+	set(MKLROOT $ENV{HOME}/intel/mkl)
+endif()
+
+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)
+option(USE_MKL "Use MKL library to build eigen" 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()
+
+include("../../external/cmake-files/QD.cmake")
+GetQD()
+
+# 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")
+
+  if (USE_MKL)
+    find_library(BLAS NAMES libblas.so libblas.dylib PATHS /usr/local/Cellar/lapack/3.9.1_1/lib /usr/lib/x86_64-linux-gnu /usr/lib/i386-linux-gnu /usr/local/Cellar/openblas/0.3.15_1/lib /usr/lib)
+    find_library(GFORTRAN NAME libgfortran.dylib PATHS /usr/local/Cellar/gcc/10.2.0_4/lib/gcc/10)
+    find_library(LAPACK NAME liblapack.dylib PATHS /usr/lib)
+    find_library(OPENMP NAME libiomp5.dylib PATHS /opt/intel/oneapi/compiler/2021.1.1/mac/compiler/lib)
+
+    include_directories (BEFORE ${MKLROOT}/include)
+    set(PROJECT_LIBS ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES} ${GFORTRAN_LIBRARIES})
+    set(MKL_LINK "-L${MKLROOT}/lib -Wl,-rpath,${MKLROOT}/lib -lmkl_intel_ilp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl")
+    add_definitions(-DEIGEN_USE_MKL_ALL)
+  else()
+    set(MKL_LINK "")
+  endif(USE_MKL)
+
+  include_directories (BEFORE ../../external)
+  include_directories (BEFORE ../../external/minimum_ellipsoid)
+  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} "-std=c++17")  # enable C++17 standard
+  set(ADDITIONAL_FLAGS "-march=native -DSIMD_LEN=0 -DTIME_KEEPING")
+  add_definitions(${CMAKE_CXX_FLAGS} "-O3 -DTIME_KEEPING"  ${ADDITIONAL_FLAGS})  # 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(volume_example volume_example.cpp)
+  target_link_libraries(volume_example QD_LIB ${MKL_LINK} ${LP_SOLVE})
+
+endif()

examples/crhmc_cooling_nonspherical_gaussians/volume_example.cpp

@@ -0,0 +1,81 @@
+// VolEsti (volume computation and sampling library)
+
+// Copyright (c) 2012-2024 Vissarion Fisikopoulos
+// Copyright (c) 2018-2024 Apostolos Chalkis
+// Copyright (c) 2024 Vladimir Necula
+
+// Contributed and/or modified by Vladimir Necula, as part of Google Summer of
+// Code 2024 program.
+
+// Licensed under GNU LGPL.3, see LICENCE file
+
+#include "generators/known_polytope_generators.h"
+#include "generators/h_polytopes_generator.h"
+#include "random_walks/random_walks.hpp"
+#include "volume/volume_cooling_nonspherical_gaussians_crhmc.hpp"
+#include "volume/volume_cooling_gaussians.hpp"
+#include <iostream>
+#include <fstream>
+#include <Eigen/Dense>
+#include <vector>
+#include "misc/misc.h"
+
+const unsigned int FIXED_SEED = 42;
+
+typedef double NT;
+typedef Cartesian<NT> Kernel;
+typedef typename Kernel::Point Point;
+typedef BoostRandomNumberGenerator<boost::mt11213b, NT, FIXED_SEED> RandomNumberGenerator;
+typedef boost::mt19937 PolyRNGType;
+typedef HPolytope<Point> HPOLYTOPE;
+
+NT nonspherical_crhmc_volume(HPOLYTOPE& polytope) {
+    int walk_len = 10;
+    NT e = 0.1;
+    RandomNumberGenerator rng;
+    // NT volume = volume_cooling_gaussians<GaussianBallWalk, RandomNumberGenerator>(polytope, e, walk_len);
+    NT volume = non_spherical_crhmc_volume_cooling_gaussians<HPOLYTOPE, RandomNumberGenerator>(polytope, rng, e, walk_len);
+    return volume;
+}
+
+NT spherical_gaussians_volume(HPOLYTOPE& polytope) {
+    int walk_len = 10;
+    NT e = 0.1;
+    RandomNumberGenerator rng;
+    NT volume = volume_cooling_gaussians<GaussianCDHRWalk, RandomNumberGenerator>(polytope, e, walk_len);
+    return volume;
+}
+
+int main() {
+
+    HPOLYTOPE cube3 = generate_cube<HPOLYTOPE>(3, false);
+    std::cout << "Cube3 \n";
+    std::cout << "Calculated Volume With Gaussian CDHR: " << spherical_gaussians_volume(cube3) << "\n";
+    std::cout << "Calculated Volume With CRHMC: " << nonspherical_crhmc_volume(cube3) << "\n";
+    std::cout << "Expected Volume: " << std::pow(2, 3) << "\n\n";
+
+    HPOLYTOPE cube4 = generate_cube<HPOLYTOPE>(4, false);
+    std::cout << "Cube4 \n";
+    std::cout << "Calculated Volume With Gaussian CDHR: " << spherical_gaussians_volume(cube4) << "\n";
+    std::cout << "Calculated Volume With CRHMC: " << nonspherical_crhmc_volume(cube4) << "\n";
+    std::cout << "Expected Volume: " << std::pow(2, 4) << "\n\n";
+
+    HPOLYTOPE skinnycube3 = generate_skinny_cube<HPOLYTOPE>(3, false);
+    std::cout << "SkinnyCube3 \n";
+    std::cout << "Calculated Volume With Gaussian CDHR: " << spherical_gaussians_volume(skinnycube3) << "\n";
+    std::cout << "Calculated Volume With CRHMC: " << nonspherical_crhmc_volume(skinnycube3) << "\n";
+    std::cout << "Expected Volume: " << 200 * std::pow(2, 2) << "\n\n";
+
+    HPOLYTOPE P3 = random_hpoly<HPOLYTOPE, PolyRNGType>(3, 12, false);
+    std::cout << "Random 3D Hpoly \n";
+    std::cout << "Calculated Volume With Gaussian CDHR: " << spherical_gaussians_volume(P3) << "\n";
+    std::cout << "Calculated Volume With CRHMC: " << nonspherical_crhmc_volume(P3) << "\n";
+    std::cout << "Expected Volume: " << "N/A" << "\n\n";
+
+    HPOLYTOPE P4 = random_hpoly<HPOLYTOPE, PolyRNGType>(4, 16, false);
+    std::cout << "Random 4D Hpoly \n";
+    std::cout << "Calculated Volume With Gaussian CDHR: " << spherical_gaussians_volume(P4) << "\n";
+    std::cout << "Calculated Volume With CRHMC: " << nonspherical_crhmc_volume(P4) << "\n";
+    std::cout << "Expected Volume: " << "N/A" << "\n\n";
+    return 0;
+}

include/ode_solvers/implicit_midpoint.hpp

@@ -172,6 +172,8 @@ struct ImplicitMidpointODESolver {
       stream << '\n';
     }
   }
+
+  NT get_eta() const { return eta; }
 };
 
 #endif

include/ode_solvers/oracle_functors.hpp

@@ -394,4 +394,77 @@ struct HessianFunctor {
 
 };
 
+struct NonSphericalGaussianFunctor {
+  template <typename NT, typename Point>
+  struct parameters {
+    Point x0;
+    NT a;
+    NT eta;
+    unsigned int order;
+    NT L;
+    NT m;
+    NT kappa;
+    Eigen::Matrix<NT, Eigen::Dynamic, Eigen::Dynamic> inv_covariance_matrix;
+
+    parameters(Point x0_, NT a_, NT eta_, Eigen::Matrix<NT, Eigen::Dynamic, Eigen::Dynamic> inv_covariance_matrix_)
+        : x0(x0_), a(a_), eta(eta_), order(2),
+          L(a_), m(a_), kappa(1),
+          inv_covariance_matrix(inv_covariance_matrix_) {};
+  };
+
+  template <typename Point>
+  struct GradientFunctor {
+    typedef typename Point::FT NT;
+    typedef std::vector<Point> pts;
+
+    parameters<NT, Point>& params;
+
+    GradientFunctor(parameters<NT, Point>& params_) : params(params_) {};
+
+    Point operator()(unsigned int const& i, pts const& xs, NT const& t) const {
+      if (i == params.order - 1) {
+          Point y = xs[0] - params.x0;
+          Eigen::Matrix<NT, Eigen::Dynamic, 1> result = -2.0 * params.a * params.inv_covariance_matrix * y.getCoefficients();
+          return Point(result);
+      } else {
+          return xs[i + 1];
+      }
+    }
+    Point operator()(Point const& x) {
+      Point y = x - params.x0;
+      Eigen::Matrix<NT, Eigen::Dynamic, 1> result = -2.0 * params.a * params.inv_covariance_matrix * y.getCoefficients();
+      return Point(result);
+    }
+  };
+
+  template <typename Point>
+  struct FunctionFunctor {
+    typedef typename Point::FT NT;
+
+    parameters<NT, Point>& params;
+
+    FunctionFunctor(parameters<NT, Point>& params_) : params(params_) {};
+
+    NT operator()(Point const& x) const {
+      Point y = x - params.x0;
+      Eigen::Matrix<NT, Eigen::Dynamic, 1> y_coeffs = y.getCoefficients();
+      return params.a * y_coeffs.dot(params.inv_covariance_matrix * y_coeffs);
+    }
+  };
+
+  template <typename Point>
+  struct HessianFunctor {
+    typedef typename Point::FT NT;
+
+    parameters<NT, Point>& params;
+
+    HessianFunctor(parameters<NT, Point>& params_) : params(params_) {};
+
+    Point operator()(Point const& x) const {
+        Eigen::Matrix<NT, Eigen::Dynamic, Eigen::Dynamic> result =   2.0 * params.a *  params.inv_covariance_matrix;
+        return Point(result);
+    }
+  };
+};
+
 #endif

include/preprocess/crhmc/opts.h

@@ -44,6 +44,7 @@ template <typename Type> class opts {
   bool DynamicWeight = true; //Enable the use of dynamic weights for each variable when sampling
   bool DynamicStepSize = true;  // Enable adaptive step size that avoids low acceptance probability
   bool DynamicRegularizer = true; //Enable the addition of a regularization term
+  bool etaInitialize = true; // enable eta initialization
   Type regularization_factor=1e-20;
   /*Dynamic step choices*/
   Type warmUpStep = 10;

include/random_walks/crhmc/additional_units/dynamic_step_size.hpp

@@ -50,7 +50,7 @@ class dynamic_step_size {
     consecutiveBadStep = IVT::Zero(simdLen);
     rejectSinceShrink = VT::Zero(simdLen);
 
-    if (options.warmUpStep > 0) {
+    if (options.warmUpStep > 0 && options.etaInitialize) {
       eta = 1e-3;
     } else {
       warmupFinished = true;

include/random_walks/crhmc/crhmc_walk.hpp

@@ -164,6 +164,9 @@ struct CRHMCWalk {
     inline void disable_adaptive(){
       update_modules=false;
     }
+    NT get_current_eta() const {
+        return solver->get_eta();
+    }
     inline void apply(RandomNumberGenerator &rng,
       int walk_length = 1,
       bool metropolis_filter = true)

include/random_walks/gaussian_helpers.hpp

@@ -9,7 +9,14 @@ NT eval_exp(Point const& p, NT const& a)
 {
     return std::exp(-a * p.squared_length());
 }
-
+template <typename Point, typename NT, typename MT>
+NT eval_exp(Point const& p, MT const& inv_covariance_matrix, NT const& a_next, NT const& a_curr)
+{
+    Eigen::Matrix<NT, Eigen::Dynamic, 1> dist_vector = p.getCoefficients();
+    NT mahalanobis_dist = dist_vector.transpose() * inv_covariance_matrix * dist_vector;
+    NT log_ratio = (a_curr - a_next) * mahalanobis_dist;
+    return std::exp(log_ratio);
+}
 
 template <typename Point, typename NT>
 NT get_max(Point const& l, Point const& u, NT const& a_i)