Parallel R with the future package¶

# Load miniconda
module load miniconda3-4.7.12.1-gcc-4.8.5-lmtvtik

# create r_parallel which installs 4
conda create -n r_parallel_tutorial r=4 r-irkernel -c conda-forge

# enables you do use conda activate
eval "$(conda shell.bash hook)"

conda activate r_parallel_tutorial

# double check you have multiple cores available

install.packages(
    c("future", "batchtools", "future.batchtools", "repurrrsive",
      "purrr","dplyr","furrr","tidyr", "snow", "ggplot2",
      "tibble", "labeling", "farver","backports"), 
    repos='http://cran.us.r-project.org', 
    Ncpus = 8 )

library(future) # needed for all examples
library(purrr) # needed for tidyvese example
library(dplyr) # needed for tidyverse example
library(furrr) # needed for tidyverse example
library(repurrrsive) # loads in data for tidyverse example
library(future.batchtools) # needed for bayesian example
library(tidyr) # needed for tidyverse example
library(snow) # needed to load Hmsc for bayesian example
library(ggplot2) # needed to plot functions
library(labeling) # needed to plot functions
library(farver) # needed to plot functions



# double-check you have access to mulitple cores 
availableCores()

# check that you can run multicore (will fail for R-studio and windows!)
supportsMulticore()

sessionInfo()

# custom function that waits a half a second and then prints current step value
slow_function <- function(step){

  # wait half a second
  Sys.sleep(.5)
    paste0("Step ", step, " completed")

}


#### running loop in serial ####
# timestamp before loop 
t1 <- proc.time()
# pre-allocate output 
output <- rep(NA, 50)

# for-loop
for (i in 1:50){
  output[i] <- slow_function(i)

}

# print the output for first 6 steps
head(output)

# timestampe after loop 
t2 <- proc.time()


print("Elapsed time (seconds) for for-loop")
# length of time
time_Elapsed <- t2[[3]] - t1[[3]]
round(time_Elapsed, 2)

#### rerunning loop in parallel using all available cores ####
# going to send iterations to each cores
# to run it only on specified number of cores, add the argument, workers. 
plan(multicore)

iterations <- 50
# rewrite loop so it works with future:
# future requires the iterative step to be saved as its own list, 
# thus we will create a list object with 50 slots for each step. 
command_set <- vector(mode = "list", length = iterations)

t1 <- proc.time()
for (i in 1:iterations){

  command_set[[i]] <- future(
    # code for each iteration within the curly brackets
    {slow_function(i)}
  ) 
}


t2 <- proc.time()


print("Elapsed time (seconds) for planning command set")
# length of time
time_Elapsed <- t2[[3]] - t1[[3]]
round(time_Elapsed, 2)

print("What one iteration looks like in the command set")
# prints first iteration
command_set[[1]]

# evaluate, distribute each iteration here 
output <- value(command_set)

# ouput is here but separates lists. 
print("list output")
head(output) # using head to limit output


# we can combine them into a vector using unlist()
print("array output")
head(unlist(output)) # using head to limit output

# plan
plan(multicore)

command_set <- vector(mode = "list", length = 50)
t1 <- proc.time()

for (i in 1:50){
  command_set[[i]] <- future(
    # but code for each iteration with curly brackets
    {slow_function(i)}
  ) 
}

# evaluate
output <- value(command_set)

# ouput is here
head(unlist(output))


# timestampe after loop 
t2 <- proc.time()


print("Elapsed time for future example:")
# length of time
time_Elapsed <- t2[[3]] - t1[[3]]
round(time_Elapsed, 2)

# create a nested dataframe with each country as a row 
# we will iterate using the map function from the purrr package
head(gap_simple)


country_nested <- 
  gap_simple %>%
  group_by(country) %>%
  nest() %>%
  ungroup() # this is necessary or furrr will be slow



# linear model we will run for each country. 
# in the real world, don't do this!! Always check model assumptions!
custom_model <- 
  function(data){
    # just to slow it down so its more obvious it's in parallel
    Sys.sleep(.2)
    lm(lifeExp ~ pop + gdpPercap + year, data = data)
  }



# Sequential form (normal purrr)
t1 <- proc.time()
model_done <- 
    country_nested %>%
    mutate(lm_obj = map(data, custom_model))

t2 <- proc.time()
print("Elapsed time for sequential purrr  example:")
# length of time
time_Elapsed <- t2[[3]] - t1[[3]]
round(time_Elapsed, 2)

# using furrr package

# use all cores available
plan(multicore) 

# this time we do no not need to the evauluation step as it is included in the future_map function.
t1 <- proc.time()
model_done <- 
    country_nested %>%
    mutate(lm_obj = future_map(data, custom_model)) # switched map to future_map!

# switch back to one core
plan(sequential)   
t2 <- proc.time()

print("Elapsed time for furrr example:")
time_Elapsed <- t2[[3]] - t1[[3]]
round(time_Elapsed, 2)

model_done$country[[13]]
model_done$lm_obj[[13]]

# bring in data for model
install.packages("coda")
install.packages("devtools") # if not yet installed
library(devtools)
install_github("hmsc-r/HMSC")
library(coda)
library(Hmsc)
load("hmsc_setup.RData")

# Setting up the model
studyDesign = data.frame(Route = XData$Route)
rL = Hmsc::HmscRandomLevel(sData=xy)
XFormula = ~ hab + poly(clim,degree = 2,raw = TRUE)
TrFormula = ~Migration + LogMass

# parameters for bayesian models. 
nChains = 4
nParallel = 4 
samples = 10 

#!/bin/bash

## Job name:
#PBS -q default
#PBS -l nodes=1:ppn=4
#PBS -l walltime=0:20:00
#PBS -N x_big_model_test_parallel_4_cores_4_chains
#PBS -j oe
#PBS -m ae


start=`date +%s`

cd $PBS_O_WORKDIR

# load R
module load r-4.0.4-gcc-10.2.0-python3-dghog6f

#Rscript -e '.libPaths("~/R/Jupyter")' -e 'batchtools::doJobCollection("<%= uri %>")'
Rscript -e 'batchtools::doJobCollection("<%= uri %>")'
end=`date +%s`

runtime=$((end-start))

echo "Runtime was $runtime seconds"

# set up that it will submit pbs scripts for each model
plan(sequential)

m = Hmsc(Y=Y, XData = XData, XFormula=XFormula, 
         phyloTree = phyloTree, TrData = TrData, 
         TrFormula = TrFormula,
         distr="probit", studyDesign=studyDesign, 
         ranLevels=list(Route=rL))
y <- list()
t1 <- proc.time()
# running 4 models, each with a different thinning value
for (thin in c(2,3,4,5)){
  y[[thin]] <- future({
  transient = 50*thin
  m = sampleMcmc(m, thin = thin, samples = samples, transient = transient,
                 nChains = nChains, initPar = "fixed effects",
                 nParallel = nParallel)

  # write model outputs to file               
  filename=file.path(paste0("Big_model_torque_chains_",as.character(nChains),"_samples_",as.character(samples),"_thin_",as.character(thin)))
  save(m,file=filename)
  }, seed = TRUE)
}
# evaluate expression
#y <- value(y) 

t2 <- proc.time()
time_Elapsed <- t2[[3]] - t1[[3]]
round(time_Elapsed, 2)

# set up that it will submit pbs scripts for each model
# calls upon the batchtools.torque.tmpl file to set the parameters for each job

plan(batchtools_torque)

m = Hmsc(Y=Y, XData = XData, XFormula=XFormula, 
         phyloTree = phyloTree, TrData = TrData, 
         TrFormula = TrFormula,
         distr="probit", studyDesign=studyDesign, 
         ranLevels=list(Route=rL))

y <- list()
t1 <- proc.time()
# running 4 models, each with a different thinning value
for (thin in c(2,3,4,5)){
  y[[thin]] <- future({
  transient = 50*thin
  m = sampleMcmc(m, thin = thin, samples = samples, transient = transient,
                 nChains = nChains, initPar = "fixed effects",
                 nParallel = nParallel)

  # write model outputs to file               
  filename=file.path(paste0("Big_model_torque_chains_",as.character(nChains),"_samples_",as.character(samples),"_thin_",as.character(thin)))
  save(m,file=filename)
  }, seed = TRUE)
}
# evaluate expression
y <- value(y) 

t2 <- proc.time()
time_Elapsed <- t2[[3]] - t1[[3]]
round(time_Elapsed, 2)

speedup_calc <- 
    function(iterations, cores, wait_time, evaluate){

    try(if(cores[1] != 1) stop("Need to run the core as 1"))

    # set up arrays for storing for-loop data
        times <- array(dim = c(length(iterations), length(cores)))
        row.names(times) <- iterations



    # running the first set outside of the future loop. 
    # skips it if since there is no set up time for serial version. 
    if (evaluate == TRUE){    

        j <- 1
        for (i in 1:length(iterations)){

            current_step <- iterations[i]
                    ### time point 1
            t1 <- proc.time()
            setup <- vector(mode = "list", length = length(iterations))

                 for (k in 1:current_step){

                        Sys.sleep(wait_time)
                    }


            ### time point
            t2 <- proc.time()

            ## calculate difference in time and convert to minutes. Write to array
            times[i, j] <- (t2[[3]] - t1[[3]]) 

                }


    }

        # loop that first run through the number of cores
        for (j in 2:length(cores)){

            plan(multicore, workers = cores[j])
            # figure out time it taakes for each iteration for j number of cores
            for (i in 1:length(iterations)){

                current_step <- iterations[i]
                ### time point 1
                t1 <- proc.time()
                setup <- vector(mode = "list", length = length(iterations))

                    for (k in 1:current_step){

                        setup[[k]] <- future({Sys.sleep(wait_time)})
                    }

                #  Allows you to run it without evaluating the code
                if (evaluate == TRUE){evaluated <-  value(setup)}

                ### time point
                t2 <- proc.time()

                ## calculate difference in time and convert to minutes. Write to array
                times[i, j] <- (t2[[3]] - t1[[3]]) 
            }


        }


        ### clean and plot all of the data
        options(repr.plot.width = 5, repr.plot.height = 2)  
        # loop is finished so now cleanaing up data
        times <- as.data.frame(times)
        colnames(times) <- cores
        times$iterations <- row.names(times)

        times$iterations <- factor(as.character(iterations), levels = iterations)
        # plot time it takes for each core/iteratiaons combo

    if (evaluate == FALSE){
        time_plot <- 
            times %>%
            pivot_longer(cols = 1:length(cores), names_to = "Cores", values_to = "Seconds") %>%
            filter(Cores != "1") %>%
            ggplot(aes(x = iterations, y = Seconds, group =  Cores, color = Cores)) + 
            geom_line() + geom_point() + 
            ggtitle("Time Needed to Set Up Task")
        print(time_plot)
        }



    if (evaluate == TRUE){   

         time_plot <- 
            times %>%
            pivot_longer(cols = 1:length(cores), names_to = "Cores", values_to = "Seconds") %>%
            ggplot(aes(x = iterations, y = Seconds, group =  Cores, color = Cores)) + 
            geom_line() + geom_point() + 
            ggtitle("Time Needed to Complete Task")
        print(time_plot)

        # speedup calcs
        speedup_df <- 1/ (times[,1:length(cores)] /times[,1] )
        speedup_df$iterations <- iterations


        speedup_plot <- 
            speedup_df %>%
            pivot_longer(cols = 1:length(cores), names_to = "Cores", values_to = "Speedup") %>%
            mutate(Cores = as.numeric(Cores)) %>%
            ggplot(aes(x = Cores, y = Speedup, group =  iterations, color = iterations)) + 
            geom_line() + geom_point() + 
            ggtitle("Speedup Calculations")

        print(speedup_plot)
        }

    }

speedup_calc(iterations = seq(from = 20, 100, by = 20), cores = c(1,2, 4, 8), wait_time = .25, evaluate = FALSE)

speedup_calc(iterations = seq(from = 20, 100, by = 20), cores = c(1,2, 4, 6, 8), wait_time = .25, evaluate = TRUE)

speedup_calc(iterations = seq(from = 100, 300, by = 100), cores = c(1,2,4,6, 8), wait_time = 1, evaluate = TRUE)

speedup_calc(iterations = seq(from = 10, 30, by = 10), cores = c(1,2,4,6, 8), wait_time = 15, evaluate = TRUE)