Fixing some R examples

man/examples/data/samples_hnr_iAB_PLT_283.txt

@@ -0,0 +1,24 @@
+-0.561720279923928 -0.980595872206212 -1.42576436471775 -1.38746518258463 -1.40328952255995 -1.32594518129933 -1.17920763919491 -1.13976931175442 -1.44704190760872 -1.80758528855703
+0.105627903577331 -0.0224760977406341 -0.033346419075399 0.0583011517226242 0.128718992901148 0.352155379240605 0.30198127014967 0.346164763337832 0.852750853220236 0.81239896064913
+0.162117647512197 0.61663482915398 0.400878181575701 0.366695795608074 0.278393392324771 0.0730857562461404 0.166434815191386 0.189330647411269 -0.125279993617079 0.126627300269534
+-0.0406942552122779 -0.576900617964391 -0.313842145078441 -0.367136023164517 -0.382100222118444 -0.354039859671709 -0.137353128426152 -0.160383552642979 -0.519155392281931 -0.475450377654029
+-0.40362251190118 -0.286473185314161 -0.0899560295412541 -0.0965770833139645 -0.172949550426655 -0.125412547613036 0.150417355582264 0.114588798098553 0.245979092609726 -0.382835955162689
+0.878246951042944 0.90514099306536 1.23213880865026 1.31832983761677 1.25986701034999 1.39995300408615 1.09368733517267 1.09091168360099 1.08437305769252 1.12958522955154
+0.498400783342133 0.587765513506474 0.253033066193859 0.370366369595461 0.449724976200575 0.365044359880781 0.387907991494453 0.383922539609155 0.188591625404994 -0.412681707489972
+0.221653648829914 0.242994434466489 -0.603741680289462 -0.481123603140553 -0.430823255733856 -0.435302829856887 -0.0456562561348847 -0.0520577661000839 -0.0385949396913763 -0.224243660239918
+-0.265084561759286 0.0918819311680734 0.33493706624183 0.3813425051495 0.403035021178183 0.477504121863839 0.532006724734535 0.545519555275899 0.766384404185368 0.831916014250039
+-0.624759045353021 -0.66962164463356 -1.13959552334022 -1.27820632084631 -1.28015141876924 -1.32693076635632 -0.909064222616852 -0.937550642541083 -0.261982584282485 -0.111118755041501
+-0.0720190591786999 -0.450621402053479 -0.863523249181356 -0.872842112502145 -1.09801753225393 -1.06373666474452 -1.11037218155514 -1.10146557900737 -0.80036781143047 -0.372168139391873
+-0.00857855953622452 0.354989948682959 1.27837817189419 1.2276458819467 1.28672546515841 1.17197787379145 1.16780775465139 1.17376822972155 1.40385625815708 1.62149324528253
+-0.229761733242758 -0.0219462023841505 -0.390687969122839 -0.302494898576672 -0.448911892337499 -0.604566355138878 -0.916877016960797 -0.906457424708054 -1.1126000897522 -1.86130665917621
+-0.589488920208699 -0.721350422835282 -0.198899309272137 -0.181407705029833 -0.118642488162701 -0.168910827232443 -0.0742514161618894 -0.0793110793398157 -0.142750010906859 0.198798486949491
+-0.0284299839046239 0.120862725181517 -0.186276678178437 -0.160980196642537 -0.293834629711641 -0.243113408459892 -0.472778364589375 -0.477458926549254 -0.153555736410121 -0.255890019644184
+-0.10591001233745 -0.362589685776349 0.325478604785015 0.378715486938594 0.26636852369326 0.161411001959766 0.60597685791033 0.610078516181616 0.572140131942936 0.511051653616906
+0.305051871362991 0.118813574717376 -0.527714492647065 -0.536667992026614 -0.580328770509275 -0.649984435188841 -0.710711338095119 -0.658502482012635 -0.902318910541349 -0.764435899578719
+0.536106973675048 1.02057897403551 0.605598476112832 0.50868937957064 0.650724511213557 0.751611789150715 0.897137728286995 0.850360925500868 1.41770846041982 1.7846819452796
+0.0564993149832159 0.463479099414294 0.460736813395207 0.474848735635176 0.426112416628774 0.390581106871229 0.793645778910023 0.806735588287548 0.844742715658839 1.09151995921733
+-0.067723281679949 0.564285422804652 0.788229511696293 0.788617644079093 0.89513355500255 0.876105020354802 0.677547485945377 0.694822145749253 1.17466687016968 0.593832946037088
+0.225364328151606 0.103238565260935 -0.0517032463701438 -0.0662367018241913 -0.0940887051913717 0.0333532293444299 0.210733515664954 0.178213696887622 -0.423217117197478 0.456888013184721
+0.0473026615725538 -0.0456028508287293 -0.177770714303978 -0.139356053949194 0.0277153774848011 -0.0567070095383378 0.00732111889985342 0.0162173446637164 0.110044897825425 -0.195869011808354
+-0.356841619301375 -0.852050102759202 -1.56144734571101 -1.61182361556797 -1.54063437363926 -1.62074713993949 -1.79194559791364 -1.8512488891935 -1.76010649453676 -1.45638353758425
+0.341165476595727 -0.498516965260338 0.0176858700381507 0.0956453734897127 0.0912542311909078 -0.0992749262782443 -0.336170653055667 -0.335003690028184 -0.639956591017614 -0.55941334662177

man/examples/generalized_hyperbolic.R

@@ -0,0 +1,66 @@
+# 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 sampling from a Generalized Hyperbolic density
+
+# Import required libraries
+library(ggplot2)
+library(volesti)
+library(numDeriv)
+library(GeneralizedHyperbolic)
+
+A = matrix(c(1, -1), ncol=1, nrow=2, byrow=TRUE)
+b = c(4,4)
+
+f <- function(x) (-log(dghyp(x)))
+grad_f <- function(x) (-ddghyp(x)/dghyp(x))
+
+x_min = matrix(0, 1, 1)
+
+# Create domain of truncation
+P <- volesti::Hpolytope(A = A, b = b)
+
+# Smoothness and strong-convexity
+L <- estimtate_lipschitz_constant(grad_f, P, 1000)
+m <- L
+
+# Warm start point from truncated Gaussian
+warm_start <- sample_points(P, n = 1, random_walk = list("nburns" = 5000), distribution = list("density" = "gaussian", "variance" = 1/L, "mode" = x_min))
+
+# Sample points
+n_samples <- 10000
+n_burns <- n_samples / 2
+
+pts <- sample_points(P, n = n_samples, random_walk = list("walk" = "HMC", "step_size" = 0.5, "nburns" = n_burns, "walk_length" = 1, "solver" = "leapfrog", "starting_point" = warm_start[,1]), distribution = list("density" = "logconcave", "negative_logprob" = f, "negative_logprob_gradient" = grad_f, "L_" = L, "m" = m))
+
+# Plot histogram
+hist(pts,
+     probability=TRUE,
+     breaks = 100,
+     border="blue",
+     main="Genrealized Hyperbolic Density with lambda = 1, alpha = 1, beta = 0, delta = 1, mu = 0",
+     xlab="Samples",
+     ylab="Density"
+)
+
+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)
+
+n_ess = min(ess(pts))
+psrf = max(psrf_univariate(pts))
+
+cat("\nEffective sample size: ", n_ess, append=TRUE)
+cat("\nPSRF: ", psrf, append=TRUE)
+cat("\n")

man/examples/metabolic.R

@@ -0,0 +1,113 @@
+# 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)
+library(R.matlab)
+
+# Sampling from logconcave density example
+
+# Helper function
+norm_vec <- function(x) sqrt(sum(x^2))
+
+
+# Load polytopes from mat file
+root <- rprojroot::find_root_file(criterion = rprojroot::has_file("DESCRIPTION"))
+metabolic_polytope_mat <- readMat(paste(root , '/man/examples/data/polytope_e_coli.mat', sep=""))
+A <- as.matrix(metabolic_polytope_mat$polytope[[1]])
+b <- as.matrix(metabolic_polytope_mat$polytope[[2]])
+center <- as.matrix(metabolic_polytope_mat$polytope[[3]])
+radius <- as.numeric(metabolic_polytope_mat$polytope[[4]])
+sigma <- 1
+dimension <- dim(A)[2]
+
+
+# Negative log-probability oracle
+f <- function(x) (norm_vec(x)^2 / (2 * sigma^2))
+
+# Negative log-probability gradient oracle
+grad_f <- function(x) (x / sigma^2)
+
+
+# Smoothness and strong-convexity
+L <- 1 / sigma^2
+m <- 1 / sigma^2
+
+# Center polytope
+b_new <- b - A %*% center
+
+# Create volesti polytope
+P <- Hpolytope(A = A, b = c(b_new))
+
+# Rounding
+#Tr <- rounding(H)
+
+#P <- Hpolytope$new(A = Tr$Mat[1:nrow(Tr$Mat), 2:ncol(Tr$Mat)], b = Tr$Mat[,1])
+
+# Center is origin (after shift)
+x_min = matrix(0, dimension, 1)
+
+# Generate samples with HNR
+start_time <- Sys.time()
+rdhr_samples <- sample_points(P, n = 10, random_walk = list("walk" = "RDHR", "nburns" = 10, "walk_length" = 1), distribution = list("density" = "gaussian", "variance" = 1/L, "mode" = x_min))
+end_time <- Sys.time()
+
+# Calculate Effective Sample size
+rdhr_ess = ess(rdhr_samples)
+min_ess <- min(rdhr_ess)
+
+# Calculate PSRF
+rdhr_psrfs = psrf_univariate(rdhr_samples)
+max_psrf = max(rdhr_psrfs)
+elapsed_time <- end_time - start_time
+
+# Print results
+cat('Min Effective Sample Size: ')
+cat(min_ess)
+cat('\n')
+cat('Maximum PSRF: ')
+cat(max_psrf)
+cat('\n')
+cat('Time per independent sample: ')
+cat(elapsed_time / min_ess)
+cat('sec')
+
+outfile <- paste(root , '/man/examples/data/samples_hnr_iAB_PLT_283.txt', sep="")
+
+write.table(rdhr_samples, file=outfile, row.names=FALSE, col.names=FALSE)
+
+start_time <- Sys.time()
+hmc_samples <- sample_points(P, n = 10, random_walk = list("walk" = "HMC", "step_size" = 0.07, "nburns" = 10, "walk_length" = 30, "solver" = "leapfrog", "starting_point" = rdhr_samples[, ncol(rdhr_samples)]), distribution = list("density" = "logconcave", "negative_logprob" = f, "negative_logprob_gradient" = grad_f, "L_" = L, "m" = m))
+end_time <- Sys.time()
+
+# Calculate Effective Sample size
+hmc_ess = ess(hmc_samples)
+min_ess <- min(hmc_ess)
+
+# Calculate PSRF
+hmc_psrfs = psrf_univariate(hmc_samples)
+max_psrf = max(hmc_psrfs)
+elapsed_time <- end_time - start_time
+
+# Print results
+cat('HMC\n')
+cat('Min Effective Sample Size: ')
+cat(min_ess)
+cat('\n')
+cat('Maximum PSRF: ')
+cat(max_psrf)
+cat('\n')
+cat('Time per independent sample: ')
+cat(elapsed_time / min_ess)
+cat('sec')
+cat('\n')

man/examples/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)
+
+# Rounding
+Tr <- rounding(H, seed = 127)
+
+P <- Hpolytope(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)
+

man/examples/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")
+}

man/examples/simple_hmc.R

@@ -0,0 +1,62 @@
+# 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)
+
+# 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(A = A, b = b)
+
+# Mode of logconcave density
+x_min <- c(-0.5)
+
+# Smoothness and strong-convexity
+L <- 2
+m <- 2
+
+# Warm start point from truncated Gaussian
+warm_start <- sample_points(P, n = 1, random_walk = list("nburns" = 5000), distribution = list("density" = "gaussian", "variance" = 1/L, "mode" = x_min))
+
+# Sample points
+n_samples <- 20000
+n_burns <- n_samples / 2
+
+pts <- sample_points(P, n = n_samples, random_walk = list("walk" = "HMC", "step_size" = 0.3, "nburns" = n_burns, "walk_length" = 3, "solver" = "leapfrog", "starting_point" = warm_start[,1]), distribution = list("density" = "logconcave", "negative_logprob" = f, "negative_logprob_gradient" = grad_f, "L_" = L, "m" = m))
+# pts <- sample_points(P, n = n_samples, random_walk = list("walk" = "HMC", "step_size" = 0.3, "nburns" = n_burns, "walk_length" = 3, "solver" = "leapfrog", "starting_point" = warm_start[,1]), distribution = list("density" = "logconcave", "mode" = x_min, "variance" = 1))
+
+# 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")