Use of SIMD computation in CRHMC walk [fixed] (#286)

R-proj/examples/logconcave/simple_crhmc.R

@@ -0,0 +1,49 @@
+# VolEsti (volume computation and sampling library)
+
+# Copyright (c) 2012-2020 Vissarion Fisikopoulos
+# Copyright (c) 2018-2020 Apostolos Chalkis
+# Copyright (c) 2020-2020 Marios Papachristou
+# Copyright (c) 2022-2022 Ioannis Iakovidis
+
+# Contributed and/or modified by Ioannis Iakovidis, as part of Google Summer of Code 2022 program.
+
+# Licensed under GNU LGPL.3, see LICENCE file
+
+# Example script for using the logconcave sampling methods
+
+# Import required libraries
+library(volesti)
+
+# Sampling from uniform density example
+
+logconcave_sample<- function(P,distribution, n_samples ,n_burns){
+  if (distribution == "uniform"){
+    f <- function(x) (0)
+    grad_f <- function(x) (0)
+    L=1
+    m=1
+    pts <- sample_points(P, n = n_samples, random_walk = list("walk" = "CRHMC", "nburns" = n_burns, "walk_length" = 1, "solver" = "implicit_midpoint"), distribution = list("density" = "logconcave", "negative_logprob" = f, "negative_logprob_gradient" = grad_f, "L_" = L, "m" = m))
+    return(max(psrf_univariate(pts, "interval")))
+  }
+  else if(distribution == "gaussian"){
+    pts <- sample_points(P, n = n_samples, random_walk = list("walk" = "CRHMC", "nburns" = n_burns, "walk_length" = 1, "solver" = "implicit_midpoint"), distribution = list("density" = "logconcave", "variance"=8))
+    return(max(psrf_univariate(pts, "interval")))
+  }
+}
+
+for (i in 1:2) {
+
+  if (i==1) {
+    distribution = 'gaussian'
+    cat("Gaussian ")
+  } else {
+    distribution = 'uniform'
+    cat("Uniform ")
+  }
+
+  P = gen_simplex(10, 'H')
+  psrf = logconcave_sample(P,distribution,5000,2000)
+  cat("psrf = ")
+  cat(psrf)
+  cat("\n")
+}

R-proj/examples/logconcave/sparse_crhmc.R

@@ -0,0 +1,98 @@
+# VolEsti (volume computation and sampling library)
+
+# Copyright (c) 2012-2020 Vissarion Fisikopoulos
+# Copyright (c) 2018-2020 Apostolos Chalkis
+# Copyright (c) 2020-2020 Marios Papachristou
+# Copyright (c) 2022-2022 Ioannis Iakovidis
+
+# Contributed and/or modified by Ioannis Iakovidis, as part of Google Summer of Code 2022 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)
+
+# Interval [-1, 1]
+A = matrix(c(1, -1), ncol=1, nrow=2, byrow=TRUE)
+b = c(2,1)
+
+# Create domain of truncation
+P <- volesti::Hpolytope$new(A, b)
+
+# Mode of logconcave density
+x_min <- c(-0.5)
+
+# Smoothness and strong-convexity
+L <- 2
+m <- 2
+
+# Sample points
+n_samples <- 80000
+n_burns <- n_samples / 2
+cat("---Sampling without hessian\n")
+pts <- sample_points(P, n = n_samples, random_walk = list("walk" = "CRHMC", "step_size" = 0.3, "nburns" = n_burns, "walk_length" = 1, "solver" = "implicit_midpoint"), distribution = list("density" = "logconcave", "negative_logprob" = f, "negative_logprob_gradient" = grad_f, "L_" = L, "m" = m))
+jpeg("histogram_without_hessian.jpg")
+# Plot histogram
+hist(pts, probability=TRUE, breaks = 100)
+
+cat("Sample mean is: ")
+sample_mean <- mean(pts)
+cat(sample_mean)
+cat("\n")
+cat("Sample variance is: ")
+sample_variance <- mean((pts - sample_mean)^2)
+cat(sample_variance)
+cat("\n")
+invisible(capture.output(dev.off()))
+
+# Negative log-probability hessian oracle
+hess_f <- function(x) (2)
+cat("---Sampling with hessian\n")
+pts <- sample_points(P, n = n_samples, random_walk = list("walk" = "CRHMC", "step_size" = 0.3, "nburns" = n_burns, "walk_length" = 1, "solver" = "implicit_midpoint"), distribution = list("density" = "logconcave", "negative_logprob" = f, "negative_logprob_gradient" = grad_f,"negative_logprob_hessian" = hess_f, "L_" = L, "m" = m))
+jpeg("histogram_with_hessian.jpg")
+# Plot histogram
+hist(pts, probability=TRUE, breaks = 100)
+
+cat("Sample mean is: ")
+sample_mean <- mean(pts)
+cat(sample_mean)
+cat("\n")
+cat("Sample variance is: ")
+sample_variance <- mean((pts - sample_mean)^2)
+cat(sample_variance)
+cat("\n")
+invisible(capture.output(dev.off()))
+
+walk="CRHMC"
+library(Matrix)
+bineq=matrix(c(10,10,10,10,10), nrow=5, ncol=1, byrow=TRUE)
+Aineq = matrix(c(1,0,-0.25,-1,2.5,1,0.4,-1,-0.9,0.5), nrow=5, ncol=2, byrow = TRUE)
+Aineq = as( Aineq, 'dgCMatrix' )
+beq=matrix(,nrow=0, ncol=1, byrow=TRUE)
+Aeq = matrix(, nrow=0, ncol=2, byrow = TRUE)
+Aeq=as( Aeq, 'dgCMatrix' )
+lb=-100000*c(1,1);
+ub=100000*c(1,1);
+cat("---Sampling the normal distribution in a pentagon\n")
+P <- volesti::sparse_constraint_problem$new(Aineq, bineq,Aeq, beq, lb, ub)
+points <- sample_points(P, n = n_samples, random_walk = list("walk" = "CRHMC", "step_size" = 0.3, "nburns" = n_burns, "walk_length" = 1, "solver" = "implicit_midpoint"), distribution = list("density" = "logconcave", "variance" = 8))
+jpeg("pentagon.jpg")
+plot(ggplot(data.frame( x=points[1,], y=points[2,] )) +
+geom_point( aes(x=x, y=y, color=walk)) + coord_fixed(xlim = c(-15,15),
+ylim = c(-15,15)) + ggtitle(sprintf("Sampling a random pentagon with walk %s", walk)))
+invisible(capture.output(dev.off()))
+write.table(points, file="pentagon.txt", row.names=FALSE, col.names=FALSE)

R-proj/man/Rcpp_sparse_constraint_problem.Rd

@@ -0,0 +1,34 @@
+\docType{class}
+\name{Rcpp_sparse_constraint_problem}
+\alias{Rcpp_sparse_constraint_problem-class}
+\alias{[,Rcpp_sparse_constraint_problem-method}
+\alias{[,Rcpp_sparse_constraint_problem,ANY,ANY,ANY-method}
+\alias{$<-,Rcpp_sparse_constraint_problem-method}
+\alias{$,Rcpp_sparse_constraint_problem-method}
+\alias{filepaths<-,Rcpp_sparse_constraint_problem-method}
+\title{
+An \code{Rcpp} class to represent sparse_constraint_problems, exposed to \code{R} via modules.
+}
+\description{
+A constraint problem is defined by a set of linear inequalities and equalities or equivalently a \eqn{d}-dimensional constraint problem is defined by a \eqn{mineq\times d} matrix Aineq and a \eqn{mineq}-dimensional vector bineq, s.t.: \eqn{Aineqx\leq bineq}, a \eqn{meq\times d} matrix Aeq and a \eqn{meq}-dimensional vector beq, s.t.: \eqn{Aeqx= beq} and two \eqn{d} vectors lb, ub such that \eqn{lb\leq x \leq ub}.}
+\details{
+\describe{
+\item{\code{Aineq} }{\eqn{mineq\times d} sparse matrix Aineq}
+
+\item{\code{bineq} }{\eqn{mineq}-dimensional vector bineq}
+
+\item{\code{Aeq} }{\eqn{meq\times d} sparse matrix Aeq}
+
+\item{\code{beq} }{\eqn{meq}-dimensional vector beq}
+
+\item{\code{lb} }{\eqn{d}-dimensional vector bineq}
+
+\item{\code{ub} }{\eqn{d}-dimensional vector bineq}
+
+\item{\code{type} }{An integer that declares the representation of the polytope. For sparse_constraint_problem the default value is 5.}
+
+\item{\code{dimension} }{The dimension of the polytope.}
+
+ }
+}
+\keyword{internal}

R-proj/man/sparse_constraint_problem.Rd

@@ -0,0 +1,26 @@
+\name{sparse_constraint_problem}
+\alias{sparse_constraint_problem}
+\title{An \code{R} class to represent sparse constraint problems.}
+
+\description{
+A constraint problem is defined by a set of linear inequalities and equalities or equivalently a \eqn{d}-dimensional constraint problem is defined by a \eqn{mineq\times d} matrix Aineq and a \eqn{mineq}-dimensional vector bineq, s.t.: \eqn{Aineqx\leq bineq}, a \eqn{meq\times d} matrix Aeq and a \eqn{meq}-dimensional vector beq, s.t.: \eqn{Aeqx\eq beq} and two \eqn{d} vectors lb, ub such that \eqn{lb\leq x \leq ub}.
+}
+\section{Fields}{
+\itemize{
+\item{\code{Aineq} }{\eqn{mineq\times d} sparse matrix Aineq}
+
+\item{\code{bineq} }{\eqn{mineq}-dimensional vector bineq}
+
+\item{\code{Aeq} }{\eqn{meq\times d} sparse matrix Aeq}
+
+\item{\code{beq} }{\eqn{meq}-dimensional vector beq}
+
+\item{\code{lb} }{\eqn{d}-dimensional vector bineq}
+
+\item{\code{ub} }{\eqn{d}-dimensional vector bineq}
+
+\item{\code{type} }{An integer that declares the representation of the polytope. For sparse_constraint_problem the default value is 5.}
+
+\item{\code{dimension} }{The dimension of the polytope.}
+
+}}

R-proj/src/Makevars

@@ -4,11 +4,16 @@ PKG_CXXFLAGS= -Wno-deprecated-declarations -lm -ldl -Wno-ignored-attributes -DBO
 
 CXX_STD = CXX11
 
-PKG_LIBS=-LRproj_externals/lp_solve -llp_solve $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)
+PKG_LIBS=-LRproj_externals/lp_solve -llp_solve -L../../external/PackedCSparse/qd -lqd $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)
 
-$(SHLIB): Rproj_externals/lp_solve/liblp_solve.a
+$(SHLIB): Rproj_externals/lp_solve/liblp_solve.a ../../external/PackedCSparse/qd/libqd.a
 
 Rproj_externals/lp_solve/liblp_solve.a:
 	@(cd Rproj_externals/lp_solve && $(MAKE) liblp_solve.a \
     CC="$(CC)" CPPFLAGS="$(CPPFLAGS)" CFLAGS="$(CFLAGS)" \
     CPICFLAGS="$(CPICFLAGS)" AR="$(AR)" RANLIB="$(RANLIB)")
+
+../../external/PackedCSparse/qd/libqd.a:
+	@(cd ../../PackedCSparse/qd && $(MAKE) libqd.a \
+		CC="$(CC)" CPPFLAGS="$(CPPFLAGS)" CFLAGS="$(CFLAGS)" \
+		CPICFLAGS="$(CPICFLAGS)" AR="$(AR)")

R-proj/src/Makevars.win

@@ -2,12 +2,17 @@ PKG_CPPFLAGS=-I../../external/boost -I../../external/LPsolve_src/run_headers -I.
 PKG_CXXFLAGS= -Wno-deprecated-declarations -lm -ldl -Wno-ignored-attributes -DBOOST_NO_AUTO_PTR -DDISABLE_NLP_ORACLES
 CXX_STD = CXX11
 
-PKG_LIBS=-LRproj_externals/lp_solve -llp_solve $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)
+PKG_LIBS=-LRproj_externals/lp_solve -llp_solve -L../../external/PackedCSparse/qd -lqd $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)
 
-$(SHLIB): Rproj_externals/lp_solve/liblp_solve.a
+$(SHLIB): Rproj_externals/lp_solve/liblp_solve.a ../../external/PackedCSparse/qd/libqd.a
 
 Rproj_externals/lp_solve/liblp_solve.a:
 	@(cd Rproj_externals/lp_solve && $(MAKE) liblp_solve.a \
     CC="$(CC)" CPPFLAGS="$(CPPFLAGS) -DUSRDLL -DINLINE=static" \
     CFLAGS="$(CFLAGS)" CPICFLAGS="$(CPICFLAGS)" AR="$(AR)" \
     RANLIB="$(RANLIB)")
+
+../../external/PackedCSparse/qd/libqd.a:
+	@(cd ../../PackedCSparse/qd && $(MAKE) libqd.a \
+		CC="$(CC)" CPPFLAGS="$(CPPFLAGS)" CFLAGS="$(CFLAGS)" \
+		CPICFLAGS="$(CPICFLAGS)" AR="$(AR)")

R-proj/src/oracle_functors_rcpp.h

@@ -17,7 +17,8 @@ enum ode_solvers {
   euler,
   runge_kutta,
   richardson,
-  collocation
+  collocation,
+  implicit_midpoint
 };
 
 // Holds Oracle Functor that wraps an R function via RCpp
@@ -91,6 +92,17 @@ struct RcppFunctor {
       }
     }
 
+    Point operator() (Point const& x) const {
+      VT y = Rcpp::as<VT>(neg_grad_f(Rcpp::wrap(x.getCoefficients())));
+
+      Point z(y);
+
+      if (negate) z = (-1.0) * z;
+
+      // Return result as Point
+      return z;
+    }
+
   };
 
   // Negative log-probability functor
@@ -118,6 +130,31 @@ struct RcppFunctor {
 
   };
 
+  // Log-probability hessian functor
+  template
+  <
+    typename Point
+  >
+  struct HessianFunctor {
+    typedef typename Point::FT NT;
+    typedef typename Point::Coeff VT;
+
+    parameters<NT> params;
+    Rcpp::Function hessian; // Negative hessian as an Rcpp::Function
+
+    HessianFunctor(
+      parameters<NT> params_,
+      Rcpp::Function hessian_) :
+      params(params_),
+      hessian(hessian_)
+    {};
+
+    Point operator() (Point const& x) const {
+      VT y= Rcpp::as<VT>(hessian(Rcpp::wrap(x.getCoefficients())));
+      return Point(y);
+    }
+
+  };
 };
 
 #endif