GeomScale · TolisChal · Jan 4, 2022 · Feb 26, 2021 · Feb 26, 2021 · Feb 26, 2021
diff --git a/R-proj/NAMESPACE b/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)
diff --git a/R-proj/R/RcppExports.R b/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.}
@@ -377,6 +376,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)

diff --git a/R-proj/src/ode_solve.cpp b/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,

diff --git a/R-proj/src/sample_points.cpp b/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,

diff --git a/R-proj/src/spectrahedron.cpp b/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;

diff --git a/R-proj/test_spectra.R b/R-proj/test_spectra.R
@@ -0,0 +1,3 @@
+library(volesti)
+
+uniform_points = sample_spectra('sdp_prob_2_6.txt', N=3000)
diff --git a/examples/spectrahedra/Boost.cmake b/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()
diff --git a/examples/spectrahedra/CMakeLists.txt b/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")

diff --git a/examples/spectrahedra/Eigen.cmake b/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()
diff --git a/examples/spectrahedra/README.md b/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. 
 

diff --git a/examples/spectrahedra/boltzmann_hmc_walk.cpp b/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;

diff --git a/examples/spectrahedra/read_write_sdpa_file.cpp b/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[]) {

diff --git a/examples/spectrahedra/solve_sdp.cpp b/examples/spectrahedra/solve_sdp.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;
 
 
 int main(int argc, char* argv[]) {

diff --git a/include/SDPAFormatManager.h b/include/SDPAFormatManager.h
@@ -244,13 +244,13 @@ class SdpaFormatManager {
     /// \param[in] is opened stream to input file
     /// \param[out] spectrahedron
     /// \param[out] objectiveFunction
-    template <typename Point>
-    void loadSDPAFormatFile(std::ifstream &is, Spectrahedron<NT, MT, VT> &spectrahedron, Point &objectiveFunction) {
+    template <typename Spectrahedron, typename Point>
+    void loadSDPAFormatFile(std::ifstream &is, Spectrahedron &spectrahedron, Point &objectiveFunction) {
         std::vector<MT> matrices;
         VT coeffs;
         loadSDPAFormatFile(is, matrices, coeffs);
         LMI<NT, MT, VT> lmi(matrices);
-        spectrahedron = Spectrahedron<NT, MT, VT>(lmi);
+        spectrahedron = Spectrahedron(lmi);
         objectiveFunction = Point(coeffs);
     }
 
@@ -260,8 +260,8 @@ class SdpaFormatManager {
     /// \param[in] is opened stream to output file
     /// \param[in] spectrahedron
     /// \param[in] objectiveFunction
-    template <typename Point>
-    void writeSDPAFormatFile(std::ostream &os, Spectrahedron<NT, MT, VT> const & spectrahedron, Point const & objectiveFunction) {
+    template <typename Spectrahedron, typename Point>
+    void writeSDPAFormatFile(std::ostream &os, Spectrahedron const & spectrahedron, Point const & objectiveFunction) {
         writeSDPAFormatFile(os, spectrahedron.getLMI().getMatrices(), objectiveFunction.getCoefficients());
     }
 };

diff --git a/include/convex_bodies/convex_body.h b/include/convex_bodies/convex_body.h
@@ -0,0 +1,107 @@
+#ifndef CONVEX_BODY_H
+#define CONVEX_BODY_H
+
+#include <limits>
+#include <iostream>
+#include <vector>
+#include <Eigen/Eigen>
+#include <functional>
+
+template <typename Point>
+class ConvexBody {
+public:
+    typedef Point                                             PointType;
+    typedef typename Point::FT                                NT;
+    typedef typename std::vector<NT>::iterator                viterator;
+    //using RowMatrixXd = Eigen::Matrix<NT, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
+    //typedef RowMatrixXd MT;
+    typedef Eigen::Matrix<NT, Eigen::Dynamic, Eigen::Dynamic> MT;
+    typedef Eigen::Matrix<NT, Eigen::Dynamic, 1>              VT;
+    typedef std::function<NT(const Point&)>                          func;
+    typedef std::function<Point(const Point&)>                       grad;
+private:
+    unsigned int        dim; //dimension
+    std::vector<func>   gs; // convex functions defining the convex body
+    std::vector<grad>   grad_gs; // convex function gradients
+    unsigned int m;
+    NT tol = NT(1e-4);
+
+
+public:
+
+    ConvexBody() : m(0) {}
+
+    ConvexBody(std::vector<func> gs_, std::vector<grad> grad_gs_, unsigned int dim_) :
+        gs(gs_), grad_gs(grad_gs_), dim(dim_)
+    {
+        m = gs.size();
+    }
+
+    unsigned int dimension() {
+        return dim;
+    }
+
+    // Compute positive line intersection (in [0, 1]) using Binary search
+    // x: starting point
+    // v: direction (ray)
+    std::pair<NT, int> line_positive_intersect(Point const& x, Point const &v) const {
+        NT t_min = NT(1);
+        int constraint = -1;
+        NT t;
+        for (unsigned int i = 0; i < m; i++) {
+            t = binary_search(x, v, gs[i]);
+            if (t < t_min) {
+                t_min = t;
+                constraint = i;
+            }
+        }
+
+        return std::make_pair(t_min, constraint);
+    }
+
+    NT binary_search(Point const &x, Point const &v, func const& f) const {
+        NT t_min = NT(0);
+        NT t_max = NT(1);
+        NT t;
+        NT value;
+
+        while (t_max - t_min > tol) {
+            t = (t_max + t_min) / 2;
+            value = f(x + t * v);
+            if (value >= -tol && value <= 0) {
+                return t;
+            } else if (value < -tol) {
+                t_min = t;
+            } else {
+                t_max = t;
+            }
+
+        }
+
+        return t;
+    }
+
+
+    // Computes unit normal at point p of the boundary
+    Point unit_normal(Point const& p, int const& constraint) const {
+        Point n = grad_gs[constraint](p);
+        return (1 / n.length()) * n;
+    }
+
+    // Computes reflection of v about point p on the boundary
+    void compute_reflection(Point &v, Point const& p, int const& constraint) const {
+        Point n = unit_normal(p, constraint);
+        v += -2 * v.dot(n) * n;
+    }
+
+    // Check if point is in K
+    int is_in(Point const& p, NT tol=NT(0)) {
+        for (func g : gs) {
+            if (g(p) > NT(-tol)) return 0;
+        }
+        return -1;
+    }
+
+};
+
+#endif
diff --git a/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h
@@ -567,7 +567,6 @@ class HPolytope {
 
 
     //First coordinate ray intersecting convex polytope
-
     std::pair<NT,NT> line_intersect_coord(Point const& r,
                                           unsigned int const& rand_coord,
                                           VT& lamdas) const
@@ -864,6 +863,8 @@ class HPolytope {
         v += -2 * v.dot(A.row(facet)) * A.row(facet);
     }
 
+    void resetFlags() {}
+
     NT log_barrier(Point &x, NT t = NT(100)) const {
       int m = num_of_hyperplanes();
       NT total = NT(0);
@@ -898,6 +899,8 @@ class HPolytope {
             v += a;
     }
 
+    void update_position_internal(NT&){}
+
     template <class bfunc, class NonLinearOracle>
     std::tuple<NT, Point, int> curve_intersect(
       NT t_prev,