Skip to content

Commit

Permalink
Initial CRHMC with Cooling Non-SphericalGaussians
Browse files Browse the repository at this point in the history
  • Loading branch information
vgnecula committed Aug 16, 2024
1 parent d424f3e commit f83f4b8
Show file tree
Hide file tree
Showing 9 changed files with 851 additions and 2 deletions.
132 changes: 132 additions & 0 deletions examples/crhmc_cooling_nonspherical_gaussians/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -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()
69 changes: 69 additions & 0 deletions examples/crhmc_cooling_nonspherical_gaussians/volume_example.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,69 @@
// 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 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";

return 0;
}
2 changes: 2 additions & 0 deletions include/ode_solvers/implicit_midpoint.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -172,6 +172,8 @@ struct ImplicitMidpointODESolver {
stream << '\n';
}
}

NT get_eta() const { return eta; }
};

#endif
73 changes: 73 additions & 0 deletions include/ode_solvers/oracle_functors.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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
1 change: 1 addition & 0 deletions include/preprocess/crhmc/opts.h
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
3 changes: 3 additions & 0 deletions include/random_walks/crhmc/crhmc_walk.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
9 changes: 8 additions & 1 deletion include/random_walks/gaussian_helpers.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
Loading

0 comments on commit f83f4b8

Please sign in to comment.