---
title: 'STAT 6984: Advanced Statistical Computing'
author:
- Jaeo Han, Byung-Jun Kim, Robert Settlage, Furong Sun and Ana (Valeria) Quevedo
date: "November 28, 2017"
output:
html_document: default
geometry: margin=2cm
setspace: doublespacing
subtitle: __GPU libraries for R__

fontsize: 12pt
---
```{r setup, eval=T, echo=F}
#misc setup stuff
graphics.off()
options(stringsAsFactors = FALSE)
knitr::opts_chunk$set(echo = T, eval=T, cache=F, tidy=T, include=FALSE, message=F, warning=F)
library.warn <- library
library <- function(package, help, pos = 2, lib.loc = NULL, character.only = FALSE,
logical.return = FALSE, warn.conflicts = FALSE, quietly = TRUE,
verbose = getOption("verbose")) {
if (!character.only) {
package <- as.character(substitute(package))
}
suppressPackageStartupMessages(library.warn(
package, help, pos, lib.loc, character.only = TRUE,
logical.return, warn.conflicts, quietly, verbose))}
options(scipen = 4, digits = 6)
# A function for captioning and referencing images
fig <- local({
i <- 0
ref <- list()
list(
cap=function(refName, text) {
i <<- i + 1
ref[[refName]] <<- i
text
},
ref=function(refName) {
ref[[refName]]
})
})
```
# Programming GPUs in R
In these exercises, you will explore the use of the GPUs as compute devices to accelerate computations using R. Our platform will be the NVidia P100 gpu's that are available on Virginia Tech's HPC cluster, NewRiver. As discussed in class, you can ssh to the cluster after getting an account with ARC (arc.vt.edu) making sure to specify you would like access to NewRiver.
## Problem 1: Interactive work with a GPU
In this problem, you will use top and nvidia-smi to profile the CPU and GPU use during two runs through a Metropolis-Hastings (MH) parallel tempering simulation.
Consider the double well potential shown in the top of Figure 1 which can be converted to a density through an exponent (bottom plots). In a normal MH simulation, you choose a start location, propose a move from that location, flip a coin to determine if you accept or reject the new location based on the probability of the coin flip vs the probability of being at the new location (+/- some indicator of how hard it was to get from the old to new location). Given some initial value, and a tunable stepsize, MH simulation will work, but will get "stuck" in one of the two hills. You could increase the stepsize, however, this can result in a "flattening" of the resulting distribution draws or non-optimum exporation of small features and in situtaions where the target has large areas of zero probability between features, still often fails to fully explore the space. Figure 2 shows MC runs exploring the bottom right distribution in Figure 1 using tuning values of 0.1 and 1. Note that we are still stuck in a single hill after widening the tuning parameter.
```{r, problem1_fig1, eval=T, echo=F, include=T, results='asis', fig.cap=fig$cap("plot1","Double well potential and unnormalized densities (gamma=1,16) density."),fig.width=5, fig.height=5}
U=function(gam,x){
gam*(x*x-1)*(x*x-1)
}
# not sure the point of the currying.
# I think it really has more to do with Darren's scala code
curried=function(gam){
#message(paste("Returning a function for gamma =",gam))
function(x) U(gam,x)
}
par(mar=c(3,2,2,1),oma=c(1,0,1,0))
layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE))
curve(U(16,x),-2,2,main="U(16)")
curve(exp(-U(1,x)),-2,2,main="exp(-U(1,x))",add = F)
curve(exp(-U(16,x)),-2,2,main="exp(-U(16,x))",add = F)
```
```{r, problem1_fig2, eval=F, echo=F, cache=F, include=T, results='asis', eval=T, fig.cap=fig$cap("plot1","Double well potential and unnormalized densities (gamma=1,16) density."),fig.width=5, fig.height=5}
set.seed(12147)
# Normal MH, adapted (sped up ;p ) from Darren Wilkinson
chains_parallel_noT_optimum=function(tune = 0.1, init = 1, temps = 16, chainLen = 1e6){
temp_len<-length(temps)
x=rep(init,temp_len)
xmat=matrix(0,chainLen,temp_len)
for (i in 1:chainLen) {
x1<-(x*x-1)*(x*x-1)
can=x+runif(temp_len,-tune,tune)
x2<-(can*can-1)*(can*can-1)
logA=temps*(x1-x2)
accept=(log(runif(temp_len))