Multiphase Monte Carlo Sampling (#147)

* implmentation of exact hmc walk

* implementation of boundary oracles for the hpolytope

* add fake boundary oracles for the other polytope representations

* implement c++ interface for exponential sampling

* improve boundary oracles and fix bugs

* add root finder for the degree two polynomial

* complete the root computation and the failure check

* initial implementation of the random walk

* improve R documentation and resolve gcc compile error

* improve R documentation

* initial implementation of hmc-leapfrog

* add the new walk to the R interface

* implement the C++ interface for the new walk

* fix compile and algorithmic errors

* improve R documentation

* improve comments and R interface

* initial implementation of the Exact HMC for gaussian dist.

* add the new walk to C++ and R interfaces

* implement boundary oracle

* improve boundary oracle for exact hmc for gaussian sampling

* remove test R script

* resolve PR reviews

* fix compile errors

* merge develop branch

* add burn in methods in exact hmc

* implement mmcs's main phase

* add window ess updater and util test functions

* delete white space in new gaussian hmc hpp file

* compare ess

* improve ess updater help functions

* resolve the reviews on the root computation

* resolve review on the header names

* improve comments and remove unused files and functions

* change a variable name in autocavariance estimator

* improve definition of macro RVOLESTI

* add an example to sample from exponential exact HMC

* add a c++ example from the gaussian exact HMC

* add two c++ examples for mmcs method, boost external files and burn in methods for aBiW

* remove RVOLESTI macro

* declare tol variable as static

* fix bug in variable declaration

* change tolerance variable declaration

* remove boost folder

* set upper bound for the number of relfections

* set upper bound for the number of relfections

* resolve review comments

* change the name of TOL to IN_INSIDE_BODY_TOLLERANCE

* implement parallel mmcs

* implement billiard walk for threads

* add example for parallel mmcs and fix compile errors

* improve comments and minor code improvements

* implement door flag for the threads in parallel mmcs

* improve parallel mmcs

* improve burnin in billiard walk

* improve the initialization of the length parameter in burnin - BiW

* improve the initialization of the length parameter in burnin - parallel BiW

* fix a bug in BiW parallel

* update parallel mmcs with openmp tools

* remove VOLESTIPY macro

* resolve PR comments

* fix c++ bugs

* resolve PR's comments

* resolve PR's comments

Co-authored-by: Vissarion Fisikopoulos <fisikop@gmail.com>

R-proj/NAMESPACE

@@ -40,4 +40,4 @@ importFrom("stats","cov")
 importFrom("utils","read.csv")
 importFrom(Rcpp,evalCpp)
 importFrom(Rcpp,loadModule)
-useDynLib(volesti, .registration=TRUE)
+useDynLib(volesti, .registration=TRUE)

R-proj/R/RcppExports.R

@@ -67,6 +67,7 @@ direct_sampling <- function(body, n, seed = NULL) {
 #' Gelman-Rubin and Brooks-Gelman Potential Scale Reduction Factor (PSRF) for each marginal
 #'
 #' @param samples A matrix that contans column-wise the sampled points from a geometric random walk.
+#' @param method A string to reauest diagnostic: (i) \code{'normal'} for psrf of Gelman-Rubin and (ii) \code{'interval'} for psrf of Brooks-Gelman.
 #'
 #' @references \cite{Gelman, A. and Rubin, D. B.,
 #' \dQuote{Inference from iterative simulation using multiple sequences,} \emph{Statistical Science,} 1992.}

R-proj/src/direct_sampling.cpp

@@ -19,7 +19,6 @@
 #include "volume/volume_sequence_of_balls.hpp"
 #include "sampling/simplex.hpp"
 
-
 //' Sample perfect uniformly distributed points from well known convex bodies: (a) the unit simplex, (b) the canonical simplex, (c) the boundary of a hypersphere or (d) the interior of a hypersphere.
 //'
 //' The \eqn{d}-dimensional unit simplex is the set of points \eqn{\vec{x}\in \R^d}, s.t.: \eqn{\sum_i x_i\leq 1}, \eqn{x_i\geq 0}. The \eqn{d}-dimensional canonical simplex is the set of points \eqn{\vec{x}\in \R^d}, s.t.: \eqn{\sum_i x_i = 1}, \eqn{x_i\geq 0}.

R-proj/src/frustum_of_simplex.cpp

@@ -3,7 +3,6 @@
 // Copyright (c) 2012-2018 Vissarion Fisikopoulos
 // Copyright (c) 2018 Apostolos Chalkis
 
-
 #include <Rcpp.h>
 #include <RcppEigen.h>
 #include "volume/exact_vols.h"
@@ -39,5 +38,4 @@ double frustum_of_simplex(Rcpp::NumericVector a, double z0){
     std::vector<double> hyp = Rcpp::as<std::vector<double> >(a);
 
     return vol_Ali(hyp, -z0, dim);
-
 }

R-proj/src/geweke.cpp

@@ -17,7 +17,6 @@
 #include <boost/random/uniform_real_distribution.hpp>
 #include "diagnostics/geweke.hpp"
 
-
 //' Geweke's MCMC diagnostic
 //'
 //' @param samples A matrix that contans column-wise the sampled points from a geometric random walk.

R-proj/src/inner_ball.cpp

@@ -1,8 +1,6 @@
-
 // Copyright (c) 2012-2018 Vissarion Fisikopoulos
 // Copyright (c) 2018 Apostolos Chalkis
 
-
 #include <Rcpp.h>
 #include <RcppEigen.h>
 #include <chrono>
@@ -97,5 +95,4 @@ Rcpp::NumericVector inner_ball(Rcpp::Reference P,
 
     vec[n] = InnerBall.second;
     return vec;
-
 }

R-proj/src/polytopes_modules.cpp

@@ -6,7 +6,6 @@
 
 #include <Rcpp.h>
 
-
 class Hpolytope {
 public:
     Hpolytope() {}
@@ -61,11 +60,9 @@ class VPinterVP {
     int type;
 };
 
-
 RCPP_MODULE(polytopes){
     using namespace Rcpp ;
 
-
     //' An exposed class to represent a H-polytope
     //'
     //' @description A H-polytope is a convex polytope defined by a set of linear inequalities or equivalently a \eqn{d}-dimensional H-polytope with \eqn{m} facets is defined by a \eqn{m\times d} matrix A and a \eqn{m}-dimensional vector b, s.t.: \eqn{Ax\leq b}.
@@ -100,7 +97,6 @@ RCPP_MODULE(polytopes){
     .field( "dimension", &Hpolytope::dimension )
     .field( "type", &Hpolytope::type );
 
-
     //' An exposed C++ class to represent a V-polytope
     //'
     //' @description A V-polytope is a convex polytope defined by the set of its vertices.
@@ -126,7 +122,6 @@ RCPP_MODULE(polytopes){
     .field( "dimension", &Vpolytope::dimension )
     .field( "type", &Vpolytope::type );
 
-
     //' An exposed C++ class to represent a zonotope
     //'
     //' @description A zonotope is a convex polytope defined by the Minkowski sum of \eqn{m} \eqn{d}-dimensional segments.
@@ -152,7 +147,6 @@ RCPP_MODULE(polytopes){
     .field( "dimension", &Zonotope::dimension )
     .field( "type", &Zonotope::type );
 
-
     //' An exposed C++ class to represent an intersection of two V-polytopes
     //'
     //' @description An intersection of two V-polytopes is defined by the intersection of the two coresponding convex hulls.

R-proj/src/rotating.cpp

@@ -99,5 +99,4 @@ Rcpp::NumericMatrix rotating (Rcpp::Reference P, Rcpp::Nullable<Rcpp::NumericMat
     MT res(TransorfMat.rows(), Rcpp::as<MT>(Mat).rows()+n);
     res << Rcpp::as<MT>(Mat).transpose(), TransorfMat;
     return Rcpp::wrap(res);
-
 }

R-proj/src/rounding.cpp

@@ -139,5 +139,4 @@ Rcpp::List rounding (Rcpp::Reference P,
     return Rcpp::List::create(Rcpp::Named("Mat") = Mat, Rcpp::Named("T") = Rcpp::wrap(std::get<0>(round_res)),
                               Rcpp::Named("shift") = Rcpp::wrap(std::get<1>(round_res)),
                               Rcpp::Named("round_value") = std::get<2>(round_res));
-
 }

R-proj/src/sample_points.cpp

@@ -24,7 +24,6 @@
 #include "ode_solvers/ode_solvers.hpp"
 #include "ode_solvers/oracle_functors_rcpp.hpp"
 
-
 enum random_walks {
   ball_walk,
   rdhr,
@@ -714,5 +713,4 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable<Rcpp::Reference> P,
     }
 
     return Rcpp::wrap(RetMat);
-
 }

R-proj/src/spectrahedron_module.cpp

@@ -20,7 +20,6 @@ class Spectrahedron {
     Spectrahedron(Rcpp::List _matrices) : matrices(_matrices) {}
 };
 
-
 RCPP_MODULE(spectrahedron){
         using namespace Rcpp ;
 

R-proj/src/volume.cpp

@@ -405,4 +405,3 @@ Rcpp::List volume (Rcpp::Reference P,
 
     return Rcpp::List::create(Rcpp::Named("log_volume") = pair_vol.first, Rcpp::Named("volume") = pair_vol.second);
 }
-

R-proj/src/zonotope_approximation.cpp

@@ -99,5 +99,4 @@ Rcpp::List zono_approx (Rcpp::Reference Z,
     }
 
     return Rcpp::List::create(Rcpp::Named("Mat") = Rcpp::wrap(Mat), Rcpp::Named("fit_ratio") = ratio);
-
 }

examples/CMakeLists.txt

@@ -1,6 +1,6 @@
 # VolEsti (volume computation and sampling library)
-# Copyright (c) 2012-2018 Vissarion Fisikopoulos
-# Copyright (c) 2018 Apostolos Chalkis
+# Copyright (c) 2012-2021 Vissarion Fisikopoulos
+# Copyright (c) 2018-2021 Apostolos Chalkis
 # Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program.
 # Licensed under GNU LGPL.3, see LICENCE file
 
@@ -112,6 +112,7 @@ else ()
   add_subdirectory(logconcave)
   #add_subdirectory(spectrahedra)
   add_subdirectory(hpolytope-volume)
+  add_subdirectory(mmcs_method)
   add_subdirectory(count-linear-extensions-using-hpolytope)
   add_subdirectory(ellipsoid-sampling)
   add_subdirectory(order-polytope-basics)

examples/mmcs_method/Boost.cmake

@@ -0,0 +1,35 @@
+function(GetBoost)
+
+	find_path(BOOST_DIR NAMES boost PATHS ../external/)
+
+	if (NOT BOOST_DIR) 
+
+    set(BOOST_URL "https://boostorg.jfrog.io/artifactory/main/release/1.76.0/source/boost_1_76_0.tar.bz2" CACHE STRING "Boost download URL")
+    set(BOOST_URL_SHA256 "f0397ba6e982c4450f27bf32a2a83292aba035b827a5623a14636ea583318c41" CACHE STRING "Boost download URL SHA256 checksum")
+
+    include(FetchContent)
+    FetchContent_Declare(
+      Boost
+      URL ${BOOST_URL}
+      URL_HASH SHA256=${BOOST_URL_SHA256}
+    )
+    FetchContent_GetProperties(Boost)
+
+    if(NOT Boost_POPULATED)
+      message(STATUS "Fetching Boost")
+      FetchContent_Populate(Boost)
+      message(STATUS "Fetching Boost - done")
+      set(BOOST_DIR ${boost_SOURCE_DIR})
+    endif()
+
+    message(STATUS "Using downloaded Boost library at ${BOOST_DIR}")
+
+  else ()
+    message(STATUS "Boost Library found: ${BOOST_DIR}")
+
+  endif()
+
+  include_directories(${BOOST_DIR})
+  include_directories(${BOOST_DIR}/boost)
+
+endfunction()

examples/mmcs_method/CMakeLists.txt

@@ -0,0 +1,124 @@
+# VolEsti (volume computation and sampling library)
+# Copyright (c) 2012-2018 Vissarion Fisikopoulos
+# Copyright (c) 2018 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 2.4.5)
+
+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("Eigen.cmake")
+GetEigen()
+
+include("Boost.cmake")
+GetBoost()
+
+# 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/Eigen)
+  include_directories (BEFORE ../../external)
+  include_directories (BEFORE ../../external/minimum_ellipsoid)
+  #include_directories (BEFORE ../../include/cartesian_geom)
+  #include_directories (BEFORE ../../include/convex_bodies)
+  include_directories (BEFORE ../../external/LPsolve_src/run_headers)
+  include_directories (BEFORE ../../external/boost)
+  #include_directories (BEFORE BOOST)
+  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)
+
+  # 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_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lgslcblas")
+  #add_definitions( "-O3 -lgsl -lm -ldl -lgslcblas" )
+
+  find_package(OpenMP REQUIRED)
+
+  add_executable (skinny_cube_10_dim skinny_cube_10_dim.cpp)
+  add_executable (random_hpoly_50_dim random_hpoly_50_dim.cpp)
+  add_executable (random_hpoly_50_dim_parallel random_hpoly_50_dim_parallel.cpp)
+
+  TARGET_LINK_LIBRARIES(skinny_cube_10_dim ${LP_SOLVE})
+  TARGET_LINK_LIBRARIES(random_hpoly_50_dim ${LP_SOLVE})
+  TARGET_LINK_LIBRARIES(random_hpoly_50_dim_parallel ${LP_SOLVE} OpenMP::OpenMP_CXX)
+
+
+endif()

examples/mmcs_method/Eigen.cmake

@@ -0,0 +1,30 @@
+function(GetEigen)
+  find_path(EIGEN_DIR NAMES Eigen PATHS ../../external)
+
+  if (NOT EIGEN_DIR) 
+    include(FetchContent)
+    FetchContent_Declare(
+      eigen
+      GIT_REPOSITORY https://gitlab.com/libeigen/eigen.git
+      GIT_TAG 3.3.9
+    )
+
+    FetchContent_GetProperties(eigen)
+
+    if(NOT eigen_POPULATED)
+      message(STATUS "Eigen library not found locally, downloading it.")
+      FetchContent_Populate(eigen)
+    endif()
+
+    set(EIGEN_DIR ${eigen_SOURCE_DIR})
+    message(STATUS "Using downloaded Eigen library at: ${EIGEN_DIR}")
+
+  else ()
+
+    message(STATUS "Eigen Library found: ${EIGEN_DIR}")
+
+  endif()
+
+  include_directories(${EIGEN_DIR})
+
+endfunction()