Support intersection between an Hpolytope and an ellipsoid

.github/workflows/R-CMD-check-ubuntu.yml

@@ -22,8 +22,8 @@ jobs:
         config:
          - {os: ubuntu-latest,  r: 'devel'}
          - {os: ubuntu-latest,  r: 'release'}
-         - {os: ubuntu-18.04,   r: 'devel'}
-         - {os: ubuntu-18.04,   r: 'release'}
+         - {os: ubuntu-20.04,   r: 'devel'}
+         - {os: ubuntu-20.04,   r: 'release'}
 
     env:
       R_REMOTES_NO_ERRORS_FROM_WARNINGS: true

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/Dirichlet/dirichlet_crhmc.R

@@ -0,0 +1,42 @@
+# VolEsti (volume computation and sampling library)
+
+# Copyright (c) 2012-2023 Vissarion Fisikopoulos
+# Copyright (c) 2018-2023 Apostolos Chalkis
+# Copyright (c) 2020-2023 Elias Tsigaridas
+
+# Licensed under GNU LGPL.3, see LICENCE file
+
+# Example script for sampling from the Dirichlet distribution with CRHMC
+
+# Import required libraries
+library(volesti)
+library(Matrix)
+
+Aeq = matrix(c(1, 1, 1), ncol=3, nrow=1, byrow=TRUE)
+Aeq = as( Aeq, 'dgCMatrix' )
+beq = matrix(c(1), ncol=1, nrow=1, byrow=TRUE)
+Aineq = matrix(, ncol=3, nrow=0, byrow=TRUE)
+Aineq = as( Aineq, 'dgCMatrix' )
+bineq = matrix(, ncol=0, nrow=0, byrow=TRUE)
+lb = c(0,0,0)
+ub = c(1, 1 ,1)
+
+# Create domain of truncation
+P <- volesti::sparse_constraint_problem$new(Aineq, bineq, Aeq, beq, lb, ub)
+
+a_vec = c(2,3,4)
+# Negative log-probability oracle
+f <- function(x) (-((a_vec-1)%*%log(x))[1])
+
+# Negative log-probability gradient oracle
+grad_f <- function(x) (-(a_vec-1)/x)
+
+n_samples <- 5000
+n_burns <- n_samples / 2
+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))

R-proj/examples/Dirichlet/dirichlet_ellipsoid_nuts.R

@@ -0,0 +1,88 @@
+# VolEsti (volume computation and sampling library)
+
+# Copyright (c) 2012-2023 Vissarion Fisikopoulos
+# Copyright (c) 2018-2023 Apostolos Chalkis
+# Copyright (c) 2020-2023 Elias Tsigaridas
+
+# Licensed under GNU LGPL.3, see LICENCE file
+
+# Example script for sampling from the Dirichlet distribution with CRHMC
+
+# Import required libraries
+library(volesti)
+library(Matrix)
+
+n = 3
+Aeq = matrix(rep(1,n), ncol=n, nrow=1, byrow=TRUE)
+beq = c(1)
+A = -diag(n)
+b = rep(0, n)
+
+# ellipsoid x'Ex <= 1
+E = matrix(c(7.8791, 3.6018, -6.9191, 3.6018, 2.2995,
+             -4.4496, -6.9191, -4.4496, 11.7993), 
+           nrow = 3, ncol = 3, byrow = TRUE)
+
+x0 = c(0.2455, 0.4287, 0.3258) #interior point
+
+# Projecting everything onto the null space of Aeq
+N_trans = pracma::nullspace(Aeq)
+N_shift = x0
+
+b = b - A %*% N_shift
+A = A %*% N_trans
+p0 = rep(0, n-1) #interior point of the projected point
+
+#the new matrix of the projected ellipsoid
+Et = t(N_trans) %*% E %*% N_trans
+#the new center of the projected ellipsoid
+center = -MASS::ginv(Et) %*% (t(N_trans) %*% E) %*% N_shift
+#the new offset of the projeced ellipsoid
+R = 1 + t(center)%*%Et%*%center - t(N_shift)%*%E%*%N_shift
+R = R[1]
+#normalizing the offset
+Et = Et / R
+
+#we shift by the ellipsoid's center so that it is given by x'Etx <= 1 
+b = b - A %*% center # shift ellipsoid's center to origin
+# interior point in the projected body
+p0 = p0 - center
+
+#now, from a point y, in the projected body, we get
+#the point oin the initial body as follows: x = N_trans%*%(y+center) + x0
+
+# Create domain of truncation
+P <- volesti::PolytopeIntersectEllipsoid$new(A, b, Et)
+
+a_vec = c(2,3,4)
+# Transformed negative log-probability oracle
+f <- function(y) (-((a_vec-1)%*%log(N_trans%*%(y+center) + x0))[1])
+
+# Transformed negative log-probability gradient oracle
+grad_f <- function(y) (-(a_vec-1)/(N_trans%*%(y+center) + x0))
+
+n_samples <- 5000
+n_burns <- n_samples / 2
+samples <- volesti::sample_points(P, n = n_samples, 
+                                  random_walk = list("walk" = "NUTS", "solver" = "leapfrog", 
+                                                     "starting_point" = p0),
+                                  distribution = list("density" = "logconcave", "negative_logprob" = f, 
+                                                      "negative_logprob_gradient" = grad_f))
+
+M = dim(samples)[2]
+# samples in the initial body following Dirichlet distribution
+random_portfolios = N_trans %*% (samples + kronecker(matrix(1, 1, M), matrix(center, ncol = 1))) + 
+  kronecker(matrix(1, 1, M), matrix(N_shift, ncol = 1))
+
+######################
+## UNIFORM SAMPLING ##
+######################
+# it's better to use Billiard walk instead of a dirichlet with ones
+
+samples <- volesti::sample_points(P, n = n_samples, 
+                                  random_walk = list("walk" = "aBiW", "starting_point" = p0))
+
+M = dim(samples)[2]
+# samples in the initial body following uniform distribution
+random_portfolios = N_trans %*% (samples + kronecker(matrix(1, 1, M), matrix(center, ncol = 1))) + 
+  kronecker(matrix(1, 1, M), matrix(N_shift, ncol = 1))

R-proj/examples/logconcave/nuts_poly_intersect_ellipsoid.R

@@ -0,0 +1,47 @@
+# VolEsti (volume computation and sampling library)
+
+# Copyright (c) 2012-2020 Vissarion Fisikopoulos
+# Copyright (c) 2018-2020 Apostolos Chalkis
+# Copyright (c) 2020-223 Elias Tsigaridas
+
+# 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)
+
+# Create domain of truncation
+A = matrix(c(1,0,0,1,-1,0,0,-1), nrow=4, ncol=2, byrow=TRUE)
+b = rep(1,4)
+E = matrix(c(0.25, 0.75, 0.75, 3.25), nrow=2, ncol=2, byrow=TRUE)
+EP = PolytopeIntersectEllipsoid$new(A, b, E)
+dim = 2
+
+x0 = rep(0, dim)
+
+# Sample points
+n_samples = 2000
+
+samples = sample_points(EP, n = n_samples, random_walk = list("walk" = "NUTS", "solver" = "leapfrog", "starting_point" = x0),
+                         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/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_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/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])