ReHMC support for spectrahedra and general convex bodies with refactored Spectrahedra classes

R-proj/NAMESPACE

@@ -18,7 +18,6 @@ export(gen_rand_vpoly)
 export(gen_rand_zonotope)
 export(gen_simplex)
 export(gen_skinny_cube)
-export(get_full_dimensional_polytope)
 export(geweke)
 export(inner_ball)
 export(loadSdpaFormatFile)
@@ -40,4 +39,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,7 +67,6 @@ 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.}
@@ -138,21 +137,6 @@ frustum_of_simplex <- function(a, z0) {
     .Call(`_volesti_frustum_of_simplex`, a, z0)
 }
 
-#' Internal rcpp function to compute the full dimensional polytope when a low dimensional is given
-#'
-#' @param P A low dimensional convex polytope in H-representation.
-#'
-#' @keywords internal
-#'
-#' @return A numerical matrix that describes the full dimensional polytope, a numerical matrix of the inverse
-#'         linear transformation that is applied on the input polytope, the numerical vector - point that the
-#'         input polytope is shifted and the product of the singular values of the matrix of the linear map 
-#'         applied on the input polytope.
-#'
-full_dimensional_polytope <- function(P) {
-    .Call(`_volesti_full_dimensional_polytope`, P)
-}
-
 #' Geweke's MCMC diagnostic
 #'
 #' @param samples A matrix that contans column-wise the sampled points from a geometric random walk.
@@ -377,6 +361,11 @@ sample_points <- function(P, n, random_walk = NULL, distribution = NULL, seed =
     .Call(`_volesti_sample_points`, P, n, random_walk, distribution, seed)
 }
 
+#' @export
+sample_spectra <- function(file = NULL, N = NULL, walk_length = NULL) {
+    .Call(`_volesti_sample_spectra`, file, N, walk_length)
+}
+
 #' Write a SDPA format file
 #'
 #' Outputs a spectrahedron (the matrices defining a linear matrix inequality) and a vector (the objective function)

R-proj/src/ode_solve.cpp

@@ -21,7 +21,7 @@
 #include "volume/volume_cooling_gaussians.hpp"
 #include "sampling/sampling.hpp"
 #include "ode_solvers/ode_solvers.hpp"
-#include "ode_solvers/oracle_functors_rcpp.hpp"
+#include "oracle_functors_rcpp.h"
 
 template <
   typename Solver,

R-proj/src/rounding.cpp

@@ -100,7 +100,7 @@ Rcpp::List rounding (Rcpp::Reference P,
             // Hpolytope
 
             if (Rcpp::as<MT>(P.field("Aeq")).rows() > 0) {
-                throw Rcpp::exception("volesti supports rounding for full dimensional polytopes. Maybe call function get_full_dimensional_polytope()");
+                throw Rcpp::exception("volesti supports rounding for full dimensional polytopes");
             } 
             Hpolytope HP(n, Rcpp::as<MT>(P.field("A")), Rcpp::as<VT>(P.field("b")));
             InnerBall = HP.ComputeInnerBall();

R-proj/src/sample_points.cpp

@@ -22,7 +22,7 @@
 #include "volume/volume_cooling_gaussians.hpp"
 #include "sampling/sampling.hpp"
 #include "ode_solvers/ode_solvers.hpp"
-#include "ode_solvers/oracle_functors_rcpp.hpp"
+#include "oracle_functors_rcpp.h"
 
 enum random_walks {
   ball_walk,

R-proj/src/spectrahedron.cpp

@@ -52,7 +52,7 @@ void writeSdpaFormatFile(Rcpp::Nullable<Rcpp::Reference> spectrahedron = R_NilVa
     typedef typename Kernel::Point Point;
     typedef boost::mt19937 RNGType;
     typedef LMI<NT, MT, VT> LMI;
-    typedef Spectrahedron<NT, MT, VT> Spectrahedron;
+    typedef Spectrahedron<Point> Spectrahedron;
 
     std::vector<MT> matrices = Rcpp::as<std::vector<MT> > (Rcpp::as<Rcpp::Reference> (spectrahedron).field("matrices"));
     LMI lmi(matrices);
@@ -91,7 +91,7 @@ Rcpp::List loadSdpaFormatFile(Rcpp::Nullable<std::string> inputFile = R_NilValue
     typedef typename Kernel::Point Point;
     typedef boost::mt19937 RNGType;
     typedef LMI<NT, MT, VT> LMI;
-    typedef Spectrahedron<NT, MT, VT> Spectrahedron;
+    typedef Spectrahedron<Point> Spectrahedron;
 
     Spectrahedron _spectrahedron;
     Point c;

examples/spectrahedra/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/spectrahedra/CMakeLists.txt

@@ -4,9 +4,30 @@
 # Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program.
 # Licensed under GNU LGPL.3, see LICENCE file
 
+cmake_minimum_required(VERSION 3.16)
+
+if(COMMAND cmake_policy)
+      cmake_policy(SET CMP0003 NEW)
+endif(COMMAND cmake_policy)
+
 
 add_executable (read_write_sdpa_file read_write_sdpa_file.cpp)
 
+include("Eigen.cmake")
+GetEigen()
+
+include("Boost.cmake")
+GetBoost()
+
+include_directories (BEFORE ../../external)
+include_directories (BEFORE ../../include/generators)
+include_directories (BEFORE ../../include/volume)
+include_directories (BEFORE ../../include)
+include_directories (BEFORE ../../include/convex_bodies)
+include_directories (BEFORE ../../include/convex_bodies/spectrahedra)
+include_directories (BEFORE ../../include/annealing)
+include_directories (BEFORE ../../include/samplers)
+include_directories (BEFORE ../../include/misc)
 
 # Find LAPACK and BLAS
 # OPENBLAS or ( ( SystemOpenblas or BLAS) and LAPACK)
@@ -26,7 +47,7 @@ IF (OPENBLAS_LIB)
         message( STATUS "ARPACK_LIB found: ${ARPACK_LIB}" )
 
         # Query gfortran to get the libgfortran path
-        FIND_LIBRARY(GFORTRAN_LIB NAMES libgfortran.so PATHS /usr/lib/gcc/x86_64-linux-gnu/8/)
+        FIND_LIBRARY(GFORTRAN_LIB NAMES libgfortran.so PATHS /usr/lib/gcc/x86_64-linux-gnu/9/)
 
         # Find libgfortran (static preferred)
         # Query gfortran to get the libgfortran path
@@ -61,7 +82,7 @@ IF (OPENBLAS_LIB)
             TARGET_LINK_LIBRARIES(boltzmann_hmc_walk ${ARPACK_LIB} ${LAPACK_LIBRARIES} ${GFORTRAN_LIB})
 
             add_executable (solve_sdp solve_sdp.cpp)
-            TARGET_LINK_LIBRARIES(solve_sdp ${ARPACK_LIB} ${LAPACK_LIBRARIES} ${GFORTRAN_LIB})
+            TARGET_LINK_LIBRARIES(solve_sdp ${GFORTRAN_LIB} ${LAPACK_LIBRARIES} ${ARPACK_LIB})
 
         ELSE()
             MESSAGE(STATUS "gfortran is required but it could not be found")

examples/spectrahedra/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.4.0
+    )
+
+    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()

examples/spectrahedra/README.md

@@ -39,7 +39,7 @@ In the folder "examples", clone this repo:
 
 ```bash
 git clone https://github.com/m-reuter/arpackpp
-cd arpackcpp
+cd arpackpp
 ```
 
 It has two scripts that should easily install the libraries:
@@ -88,7 +88,7 @@ The contents of the file are:
  1 -0 -0
 ```
 
- 
+
 It represents a spectrahedron in 2 dimensions, described by a linear matrix inequality with
 3x3 matrices, and a linear objective function. 
 

examples/spectrahedra/boltzmann_hmc_walk.cpp

@@ -28,7 +28,7 @@ typedef Eigen::Matrix<NT, Eigen::Dynamic, 1> VT;
 typedef Eigen::Matrix <NT, Eigen::Dynamic, Eigen::Dynamic> MT;
 typedef Cartesian <NT> Kernel;
 typedef typename Kernel::Point Point;
-typedef Spectrahedron <NT, MT, VT> SPECTRAHEDRON;
+typedef Spectrahedron <Point> SPECTRAHEDRON;
 typedef BoostRandomNumberGenerator<boost::mt19937, NT> RNGType;
 typedef BoltzmannHMCWalk::Walk<SPECTRAHEDRON, RNGType> HmcWalk;
 typedef BoltzmannHMCWalk::Walk<SPECTRAHEDRON, RNGType>::Settings HmcWalkSettings;

examples/spectrahedra/read_write_sdpa_file.cpp

@@ -26,7 +26,7 @@ typedef Eigen::Matrix<NT, Eigen::Dynamic, 1> VT;
 typedef Eigen::Matrix <NT, Eigen::Dynamic, Eigen::Dynamic> MT;
 typedef Cartesian <NT> Kernel;
 typedef typename Kernel::Point Point;
-typedef Spectrahedron <NT, MT, VT> SPECTRAHEDRON;
+typedef Spectrahedron <Point> SPECTRAHEDRON;
 
 
 int main(int argc, char* argv[]) {