resolve conflicts from develop

R-proj/R/RcppExports.R

@@ -361,11 +361,6 @@ 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/examples/logconcave/nuts_rand_poly.R

@@ -0,0 +1,55 @@
+# VolEsti (volume computation and sampling library)
+
+# Copyright (c) 2012-2020 Vissarion Fisikopoulos
+# Copyright (c) 2018-2020 Apostolos Chalkis
+# Copyright (c) 2020-2020 Marios Papachristou
+
+# Contributed and/or modified by Marios Papachristou, as part of Google Summer of Code 2020 program.
+
+# Licensed under GNU LGPL.3, see LICENCE file
+
+# Example script for using the logconcave sampling methods
+
+# Import required libraries
+library(ggplot2)
+library(volesti)
+
+# Sampling from logconcave density example
+
+# Helper function
+norm_vec <- function(x) sqrt(sum(x^2))
+
+# Negative log-probability oracle
+f <- function(x) (norm_vec(x)^2 + sum(x))
+
+# Negative log-probability gradient oracle
+grad_f <- function(x) (2 * x + 1)
+
+dimension <- 50
+facets <- 200
+
+# Create domain of truncation
+H <- gen_rand_hpoly(dimension, facets, seed = 15)
+
+# Rounding
+Tr <- rounding(H, seed = 127)
+
+P <- Hpolytope$new(A = Tr$Mat[1:nrow(Tr$Mat), 2:ncol(Tr$Mat)], b = Tr$Mat[,1])
+
+x_min = matrix(0, dimension, 1)
+
+# Warm start point from truncated Gaussian
+warm_start <- sample_points(P, n = 1, random_walk = list("nburns" = 5000), distribution = list("density" = "gaussian", "variance" = 1/2, "mode" = x_min))
+
+# Sample points
+n_samples <- 20000
+
+samples <- sample_points(P, n = n_samples, random_walk = list("walk" = "NUTS", "solver" = "leapfrog", "starting_point" = warm_start[,1]),
+                         distribution = list("density" = "logconcave", "negative_logprob" = f, "negative_logprob_gradient" = grad_f))
+
+# Plot histogram
+hist(samples[1,], probability=TRUE, breaks = 100)
+
+psrfs <- psrf_univariate(samples)
+n_ess <- ess(samples)
+

R-proj/examples/logconcave/simple_hmc_rand_poly.R

@@ -29,10 +29,10 @@ dimension <- 50
 facets <- 200
 
 # Create domain of truncation
-H <- gen_rand_hpoly(dimension, facets)
+H <- gen_rand_hpoly(dimension, facets, seed = 15)
 
 # Rounding
-Tr <- rounding(H)
+Tr <- rounding(H, seed = 127)
 
 P <- Hpolytope$new(A = Tr$Mat[1:nrow(Tr$Mat), 2:ncol(Tr$Mat)], b = Tr$Mat[,1])
 

R-proj/src/sample_points.cpp

@@ -36,6 +36,7 @@ enum random_walks {
   brdhr,
   bcdhr,
   hmc,
+  nuts,
   gaussian_hmc,
   exponential_hmc,
   uld,
@@ -219,6 +220,26 @@ void sample_from_polytope(Polytope &P, int type, RNGType &rng, PointList &randPo
       execute_crhmc<Polytope, RNGType, PointList, NegativeGradientFunctor,NegativeLogprobFunctor, HessianFunctor, CRHMCWalk, 8>(P, rng, randPoints, walkL, numpoints, nburns, F, f, h);
       break;
       }
+    case nuts:
+
+      logconcave_sampling <
+              PointList,
+              Polytope,
+              RNGType,
+              NutsHamiltonianMonteCarloWalk,
+              NT,
+              Point,
+              NegativeGradientFunctor,
+              NegativeLogprobFunctor,
+              LeapfrogODESolver <
+                Point,
+                NT,
+                Polytope,
+                NegativeGradientFunctor
+              >
+            >(randPoints, P, rng, walkL, numpoints, StartingPoint, nburns, *F, *f);
+
+        break;
     case uld:
 
       logconcave_sampling <
@@ -250,7 +271,17 @@ void sample_from_polytope(Polytope &P, int type, RNGType &rng, PointList &randPo
 //' @param n The number of points that the function is going to sample from the convex polytope.
 //' @param random_walk Optional. A list that declares the random walk and some related parameters as follows:
 //' \itemize{
-//' \item{\code{walk} }{ A string to declare the random walk: i) \code{'CDHR'} for Coordinate Directions Hit-and-Run, ii) \code{'RDHR'} for Random Directions Hit-and-Run, iii) \code{'BaW'} for Ball Walk, iv) \code{'BiW'} for Billiard walk, v) \code{'dikin'} for dikin walk, vi) \code{'vaidya'} for vaidya walk, vii) \code{'john'} for john walk, viii) \code{'BCDHR'} boundary sampling by keeping the extreme points of CDHR or ix) \code{'BRDHR'} boundary sampling by keeping the extreme points of RDHR x) \code{'HMC'} for Hamiltonian Monte Carlo (logconcave densities) xi) \code{'ULD'} for Underdamped Langevin Dynamics using the Randomized Midpoint Method xii) \code{'ExactHMC'} for exact Hamiltonian Monte Carlo with reflections (spherical Gaussian or exponential distribution). The default walk is \code{'aBiW'} for the uniform distribution or \code{'CDHR'} for the Gaussian distribution and H-polytopes and \code{'BiW'} or \code{'RDHR'} for the same distributions and V-polytopes and zonotopes.}
+//' \item{\code{walk} }{ A string to declare the random walk: i) \code{'CDHR'} for Coordinate Directions Hit-and-Run, 
+//' ii) \code{'RDHR'} for Random Directions Hit-and-Run, iii) \code{'BaW'} for Ball Walk, iv) \code{'BiW'} for Billiard walk, 
+//' v) \code{'dikin'} for dikin walk, vi) \code{'vaidya'} for vaidya walk, vii) \code{'john'} for john walk, 
+//' viii) \code{'BCDHR'} boundary sampling by keeping the extreme points of CDHR or ix) \code{'BRDHR'} boundary sampling by keeping the extreme points of RDHR, 
+//' x) \code{'NUTS'} for NUTS Hamiltonian Monte Carlo sampler (logconcave densities), xi) \code{'HMC'} for Hamiltonian Monte Carlo  (logconcave densities), 
+//' xii) CRHMC for Riemannian HMC with H-polytope constraints (uniform and general logconcave densities), 
+//' xiii) \code{'ULD'} for Underdamped Langevin Dynamics using the Randomized Midpoint Method (logconcave densities), 
+//' xiii) \code{'ExactHMC'} for exact Hamiltonian Monte Carlo with reflections (spherical Gaussian or exponential distribution). 
+//' The default walk is \code{'aBiW'} for the uniform distribution, \code{'CDHR'} for the Gaussian distribution and H-polytopes and 
+//' \code{'BiW'} or \code{'RDHR'} for the same distributions and V-polytopes and zonotopes. \code{'NUTS'} is the default sampler for logconcave densities and \code{'CRHMC'} 
+//' for logconcave densities with H-polytope and sparse constrainted problems.}
 //' \item{\code{walk_length} }{ The number of the steps per generated point for the random walk. The default value is \eqn{1}.}
 //' \item{\code{nburns} }{ The number of points to burn before start sampling. The default value is \eqn{1}.}
 //' \item{\code{starting_point} }{ A \eqn{d}-dimensional numerical vector that declares a starting point in the interior of the polytope for the random walk. The default choice is the center of the ball as that one computed by the function \code{inner_ball()}.}
@@ -421,18 +452,24 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable<Rcpp::Reference> P,
         Rcpp::Function negative_logprob = Rcpp::as<Rcpp::List>(distribution)["negative_logprob"];
         Rcpp::Function negative_logprob_gradient = Rcpp::as<Rcpp::List>(distribution)["negative_logprob_gradient"];
 
-        NT L_, m;
+        NT L_ = 1, m = 1;
 
         if (Rcpp::as<Rcpp::List>(distribution).containsElementNamed("L_")) {
             L_ = Rcpp::as<NT>(Rcpp::as<Rcpp::List>(distribution)["L_"]);
+            if (L_ <= NT(0)) {
+                throw Rcpp::exception("The smoothness constant must be positive");
+            }
         } else {
-            throw Rcpp::exception("The smoothness constant is absent");
+            L_ = -1;
         }
 
         if (Rcpp::as<Rcpp::List>(distribution).containsElementNamed("m")) {
             m = Rcpp::as<NT>(Rcpp::as<Rcpp::List>(distribution)["m"]);
+            if (m <= NT(0)) {
+                throw Rcpp::exception("The strong-convexity constant must be positive");
+            }
         } else {
-            throw Rcpp::exception("The strong-convexity constant is absent");
+            m = -1;
         }
 
         if (Rcpp::as<Rcpp::List>(random_walk).containsElementNamed("step_size")) {
@@ -456,7 +493,6 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable<Rcpp::Reference> P,
             throw Rcpp::exception("Invalid ODE solver specified. Aborting.");
           }
          } else {
-          Rcpp::warning("Solver set to leapfrog.");
           solver = leapfrog;
         }
         // Create functors
@@ -511,18 +547,12 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable<Rcpp::Reference> P,
                 throw Rcpp::exception("Exponential sampling is supported only for H-polytopes");
             }
             walk = exponential_hmc;
+        } else if (logconcave) {
+            walk = (type == 5) ? crhmc : nuts;
         } else if (gaussian) {
-            if (type == 1) {
-                walk = cdhr;
-            } else {
-                walk = rdhr;
-            }
+            walk = (type == 1) ? cdhr : rdhr;
         } else {
-            if (type == 1) {
-                walk = accelarated_billiard;
-            } else {
-                walk = billiard;
-            }
+            walk = (type == 1) ? accelarated_billiard : billiard;
         }
     } else if (Rcpp::as<std::string>(Rcpp::as<Rcpp::List>(random_walk)["walk"]).compare(std::string("CDHR")) == 0) {
         walk = cdhr;
@@ -587,7 +617,12 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable<Rcpp::Reference> P,
         }
     } else if (Rcpp::as<std::string>(Rcpp::as<Rcpp::List>(random_walk)["walk"]).compare(std::string("HMC")) == 0) {
         if (!logconcave) throw Rcpp::exception("HMC is not supported for non first-order sampling");
+        if (F->params.L < 0) throw Rcpp::exception("The smoothness constant is absent");
+        if (F->params.m < 0) throw Rcpp::exception("The strong-convexity constant is absent");
         walk = hmc;
+    } else if (Rcpp::as<std::string>(Rcpp::as<Rcpp::List>(random_walk)["walk"]).compare(std::string("NUTS")) == 0) {
+        if (!logconcave) throw Rcpp::exception("NUTS is not supported for non first-order sampling");
+        walk = nuts;
     } else if (Rcpp::as<std::string>(Rcpp::as<Rcpp::List>(random_walk)["walk"]).compare(std::string("ULD")) == 0) {
         if (!logconcave) throw Rcpp::exception("ULD is not supported for non first-order sampling");
         walk = uld;
@@ -747,7 +782,7 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable<Rcpp::Reference> P,
     if (numpoints % 2 == 1 && (walk == brdhr || walk == bcdhr)) numpoints--;
     MT RetMat(dim, numpoints);
     unsigned int jj = 0;
-
+    
     for (typename std::list<Point>::iterator rpit = randPoints.begin(); rpit!=randPoints.end(); rpit++, jj++) {
         if (gaussian) {
             RetMat.col(jj) = (*rpit).getCoefficients() + mode.getCoefficients();

examples/logconcave/CMakeLists.txt

@@ -20,7 +20,7 @@ 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" ON)
+option(USE_MKL "Use MKL library to build eigen" OFF)
 
 if(DISABLE_NLP_ORACLES)
   add_definitions(-DDISABLE_NLP_ORACLES)
@@ -73,16 +73,24 @@ else ()
     include_directories(BEFORE /usr/include/eigen3)
 endif(BUILTIN_EIGEN)
 
+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)
+
 if (USE_MKL)
-  find_library(BLAS NAME libblas.so PATHS /usr/lib/x86_64-linux-gnu /usr/lib/i386-linux-gnu)
+  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})
+  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)
 
 # Find lpsolve library
-find_library(LP_SOLVE NAMES liblpsolve55.so PATHS /usr/lib/lp_solve)
+find_library(LP_SOLVE NAMES liblpsolve55.so liblpsolve55.dylib PATHS /usr/lib/lp_solve /usr/local/lib)
 
 if (NOT LP_SOLVE)
   message(FATAL_ERROR "This program requires the lp_solve library, and will not be compiled.")
@@ -145,6 +153,6 @@ else ()
   TARGET_LINK_LIBRARIES(simple_hmc ${LP_SOLVE} ${BLAS} ${MKL_LINK} )
   TARGET_LINK_LIBRARIES(exponential_exact_hmc ${LP_SOLVE} ${BLAS} ${MKL_LINK}  )
   TARGET_LINK_LIBRARIES(gaussian_exact_hmc ${LP_SOLVE} ${BLAS} ${MKL_LINK}  )
-  
+
 
 endif()

external/cmake-files/LPSolve.cmake

@@ -26,7 +26,7 @@ function(GetLPSolve)
       message(STATUS "lp_solve library found: ${LP_SOLVE_DIR}")
 
   endif()
-
+  
   include_directories(${LP_SOLVE_DIR})
 
 endfunction()

include/cartesian_geom/point.h

@@ -81,6 +81,7 @@ class point
     void set_dimension(const unsigned int dim)
     {
         d = dim;
+        coeffs.setZero(d);
     }
 
     void set_coord(const unsigned int i, FT coord)