View source for R parallel

= Introduction =

* [http://shop.oreilly.com/product/0636920021421.do Example code] for the book Parallel R by McCallum and Weston.
* [http://www.win-vector.com/blog/2016/01/parallel-computing-in-r/ A gentle introduction to parallel computing in R]
* [http://www.stat.berkeley.edu/scf/paciorek-distribComp.pdf An introduction to distributed memory parallelism in R and C]
* [http://danielmarcelino.com/parallel-processing/Parallel Processing: When does it worth?]
* [https://arxiv.org/pdf/1912.11144.pdf Parallel Computing With R: A Brief Review] by Dirk Eddelbuettel (updated 2020).

= Security warning from Windows/Mac =
It seems it is safe to choose 'Cancel' when Windows Firewall tried to block R program when we use '''makeCluster()''' to create a socket cluster.
<syntaxhighlight lang='rsplus'>
library(parallel)
cl <- makeCluster(2)
clusterApply(cl, 1:2, get("+"), 3)
stopCluster(cl)
</syntaxhighlight>
[[File:WindowsSecurityAlert.png|100px]]  [[File:RegisterDoParallel mac.png|150px]]

If we like to see current firewall settings, just click Windows Start button, search 'Firewall' and choose 'Windows Firewall with Advanced Security'. In the 'Inbound Rules', we can see what programs (like, R for Windows GUI front-end, or Rserve) are among the rules. These rules are called 'private' in the 'Profile' column. Note that each of them may appear twice because one is 'TCP' protocol and the other one has a 'UDP' protocol.

= Detect number of cores and memory consideration =
<syntaxhighlight lang='rsplus'>
parallel::detectCores()  # not good
future::availableCores() # good
</syntaxhighlight>
Don't use the default option getOption("mc.cores", 2L) (PS it only returns 2.) in mclapply() unless you are a developer for a package. 

However, it is a different story when we run the R code in HPC cluster. Read the discussion [https://stackoverflow.com/questions/28954991/whether-to-use-the-detectcores-function-in-r-to-specify-the-number-of-cores-for Whether to use the detectCores function in R to specify the number of cores for parallel processing?]

On NIH's biowulf, even I specify an interactive session with 4 cores, the parallel::detectCores() function returns 56. This number is the same as the output from the bash command '''grep processor /proc/cpuinfo''' or (better) '''lscpu'''. The '''free -hm''' also returns a full 125GB size instead of my requested size (4GB by default). The '''future::availableCores()''' fixes the problem. 

See [https://hpc.nih.gov/apps/R.html#parallel Biowulf's R webpage] for a detailed instruction and my note [[Biowulf#Running_R_scripts|Running R scripts]] for selecting the number of cores on biowulf/non-biowulf (different OS).

== Identify OS ==
See [[R#Distinguish_Windows_and_Linux.2FMac|Distinguish Windows, Linux, Mac]]

== Memory ==
* [https://stackoverflow.com/a/57035805 Alternative to R's `memory.size()` in linux?] pryr::mem_used() will return the memory currently used by R.
* '''system('grep MemTotal /proc/meminfo', intern = TRUE)''' will show the total physical RAM on Linux OS (mac does not work).
* [https://stackoverflow.com/q/46248139 Monitor memory usage with R parallel]. [https://rstudio.github.io/profvis/ profvis] package from RStudio.
** It works on RStudio or a regular R session
** The output is an HTML file
** We can change the interval period. Its default is 0.01 seconds or 100 ms. '''profvis(expr = NULL, interval = 0.01, ...)'''
** When I test profvis() on a function which calls clusterExport() and parSapply(), I can see the time spent on each of them. My example of running coxph on a set of genes.
** Question: why the sum of time from the HTML output not close at all to system.time() when parLapply() is used?
*** [https://github.com/rstudio/profvis/issues/70 Running time under-estimated with parallel::parLapply()]. ''The profiler only considers the main thread. Which means that any thread/process spawned from the main rsession process, will not be accounted for in the execution time. Only the thread that spawns the workers are timed''
*** [https://github.com/rstudio/profvis/issues/72 Profile non-cpu time?]
** [https://adv-r.hadley.nz/perf-measure.html Advanced R] > Measuring performance
** [https://bookdown.org/csgillespie/efficientR/performance.html#performance-profvis Efficient R programming] > Code profiling.
* Rprof() and [https://www.rdocumentation.org/packages/utils/versions/3.6.2/topics/summaryRprof summaryRprof()]
:<syntaxhighlight lang='rsplus'>
Rprof(tf <- "rprof.txt", memory.profiling = TRUE)
# the code you want to profile must be in between
Rprof (NULL) ; print(summaryRprof(tf))
</syntaxhighlight>

= Socket or Fork =
[https://statcompute.wordpress.com/2019/06/30/parallel-r-socket-or-fork/ Parallel R: Socket or Fork]
* For the '''fork''', each parallel thread is a complete duplication of the master process with the shared environment, including objects or variables defined prior to the kickoff of parallel threads. Therefore, it runs '''fast'''. However, the major limitation is that the fork doesn’t work on the Windows system.
* The '''socket''' works on all operating systems. Each thread runs separately without sharing objects or variables, which can only be passed from the master process explicitly. As a result, it runs slower due to the communication overhead.
<syntaxhighlight lang='rsplus'>
install.packages(c("rbenchmark", "nycflights13"))
data(flights, package = 'nycflights13')
dim(flights) # [1] 336776     19
ex <- expression(carrier == "UA" & origin == "EWR" & day == 1 & is.na(arr_time))

parFilter <- function(df, ex, type) {
  cn <- parallel::detectCores() - 1
  cl <- parallel::makeCluster(cn, type = type)
  ### DIVIDE THE DATAFRAME BASED ON # OF CORES
  sp <- parallel::parLapply(cl, parallel::clusterSplit(cl, seq(nrow(df))),
                            function(c_) df[c_,])
  ### PASS THE OBJECT FROM MASTER PROCESS TO EACH NODE
  parallel::clusterExport(cl, "ex")
  ### EXTRACT ROW INDEX ON EACH NODE
  id <- Reduce(c, parallel::parLapply(cl, sp,
                                      function(s_) with(s_, eval(ex))))
  parallel::stopCluster(cl)
  return(df[which(id),])
}

rbenchmark::benchmark(replications = 10, order = "elapsed", relative = "elapsed",
                      columns = c("test", "replications", "elapsed", "relative"),
                      "  FORK" = parFilter(flights, ex, "FORK"),
                      "SOCKET" = parFilter(flights, ex, "PSOCK")
)  # Phenom(tm) II X6 1055T
    test replications elapsed relative
1   FORK           10  42.154    1.000
2 SOCKET           10  68.901    1.635
</syntaxhighlight>

== makeCluster() for socket connection in R 4.0 ==
The function hangs on R 4.0.0 + RStudio. See
* https://github.com/rstudio/rstudio/issues/6692 for a simple solution '''cl <- makeCluster(ncores, setup_strategy = "sequential") ''' 
* The article [https://developer.r-project.org/Blog/public/2020/03/17/socket-connections-update/index.html Socket Connections Update] for R 4.0

= How can I print when using %dopar% =
https://stackoverflow.com/a/15078540

Create and register your cluster with something like:
<syntaxhighlight lang='rsplus'>
library(doParallel)
cl<-makeCluster(4, outfile = "debugParallel.txt")
registerDoParallel(cl)
foreach(x=list(1, 2, "a"))  %dopar%  
{
  print(x)
}
stopCluster(cl)
</syntaxhighlight>
And the ''debugParallel.txt'' will look like
<pre>
starting worker pid=17718 on localhost:11704 at 15:10:29.665
starting worker pid=17733 on localhost:11704 at 15:10:29.850
starting worker pid=17747 on localhost:11704 at 15:10:30.025
starting worker pid=17761 on localhost:11704 at 15:10:30.183
[1] 1
[1] 2
[1] "a"
</pre>
Setting outfile to the empty string ("") will prevent doParallel (or snow whatever) from redirecting the output, often resulting in the output from your print messages showing up on the terminal of the master process. Your foreach loop doesn't need to change at all.

= parallel package (including parLapply, parSapply) =
Parallel package was included in R 2.14.0. It is derived from the snow and multicore packages and provides many of the same functions as those packages.

The parallel package provides several *apply functions for R users to quickly modify their code using parallel computing.

* makeCluster(makePSOCKcluster, makeForkCluster), stopCluster. Other cluster types are passed to package '''snow'''.
* '''clusterCall''', clusterEvalQ: source R files and/or load libraries
* clusterSplit
* '''clusterApply''', '''clusterApplyLB''' (vs the foreach package)
* '''clusterExport''': export variables. We may need to apply the option '''envir=environment()''' if [https://stackoverflow.com/a/12023454 the objects to be exported are local variables].
* clusterMap
* [https://stat.ethz.ch/R-manual/R-devel/library/parallel/html/mclapply.html parallel::mclapply()] and [https://stat.ethz.ch/R-manual/R-devel/library/parallel/html/clusterApply.html parallel::parLapply()]
* parLapply, parSapply, parApply, parRapply, parCapply. Note that 
** '''parSapply()''' can be used to as a parallel version of the replicate() function. See [https://stackoverflow.com/questions/19281010/simplest-way-to-do-parallel-replicate?answertab=active#tab-top this example] or [https://jozef.io/r920-christmas-praise-2019/#base-package-parallel this example]. 
** An iteration parameter needs to be added to the first parameter of the main function.
* parLapplyLB, parSapplyLB (load balance version)
* clusterSetRNGStream, nextRNGStream, nextRNGSubStream

Examples (See ?[http://www.inside-r.org/r-doc/parallel/clusterApply clusterApply])
<syntaxhighlight lang='rsplus'>
library(parallel)
cl <- makeCluster(2, type = "SOCK")
clusterApply(cl, 1:2, function(x) x*3)    # OR clusterApply(cl, 1:2, get("*"), 3)
# [[1]]
# [1] 3
#
# [[2]]
# [1] 6
parSapply(cl, 1:20, get("+"), 3)
#  [1]  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
stopCluster(cl)
</syntaxhighlight>
An example of using clusterCall() or clusterEvalQ()
<syntaxhighlight lang='rsplus'>
library(parallel)

cl <- makeCluster(4)
clusterCall(cl, function() { 
  source("test.R") 
})
# clusterEvalQ(cl, {
#    source("test.R")
# })

## do some parallel work
stopCluster(cl)
</syntaxhighlight>

== mclapply() (Windows OS only supports 1 core) vs parLapply()/parSapply() (all OSs) ==
* Be providing the number of cores in mclapply() using '''mc.cores''' argument (2 is used by default)
* set.seed() can not fix seed in mclapply(); we need to use set.seed(, "L'Ecuyer-CMRG"). But be careful on the side-effect of using '''L'Ecuyer-CMRG''' seed. 
* Chapter 4 of [https://www.oreilly.com/library/view/parallel-r/9781449317850/ch04.html Parallel R] by Stephen Weston, Q. Ethan McCallum
* [https://irudnyts.github.io//setting-a-seed-in-r-when-using-parallel-simulation/ Setting a seed in R, when using parallel simulation]
* '''[https://stackoverflow.com/questions/15070377/r-doesnt-reset-the-seed-when-lecuyer-cmrg-rng-is-used R doesn't reset the seed when “L'Ecuyer-CMRG” RNG is used?]''' <syntaxhighlight lang='rsplus'>
library(parallel)
system.time(mclapply(1:1e4L, function(x) rnorm(x)))
system.time(mclapply(1:1e4L, function(x) rnorm(x), mc.cores = 4))

set.seed(1234)
mclapply(1:3, function(x) rnorm(x))
set.seed(1234)
mclapply(1:3, function(x) rnorm(x)) # cannot reproduce the result

mclapply(1:3, function(x) {set.seed(x); rnorm(x)}) # Can reproduce the result

set.seed(123, "L'Ecuyer")
mclapply(1:3, function(x) rnorm(x))
mclapply(1:3, function(x) rnorm(x)) # results are not changed once we have run set.seed( , "L'Ecuyer")

set.seed(1234)                      # use set.seed() in order to get a new reproducible result
mclapply(1:3, function(x) rnorm(x)) 
mclapply(1:3, function(x) rnorm(x)) # results are not changed

# an example from ?mclapply
library(boot)
cd4.rg <- function(data, mle) MASS::mvrnorm(nrow(data), mle$m, mle$v)
cd4.mle <- list(m = colMeans(cd4), v = var(cd4))
mc <- getOption("mc.cores", 2)
run1 <- function(...) boot(cd4, corr, R = 500, sim = "parametric",
                           ran.gen = cd4.rg, mle = cd4.mle)
## To make this reproducible:
set.seed(123, "L'Ecuyer")
res <- mclapply(seq_len(mc), run1)
cd4.boot <- do.call(c, res)
boot.ci(cd4.boot,  type = c("norm", "basic", "perc"),
        conf = 0.9, h = atanh, hinv = tanh)
</syntaxhighlight>
* An example of applying to the permutation test (LS test for GSEA)
:<syntaxhighlight lang='rsplus'>
ngenes <- 10000
nresamp <- 1e5
genesetSize <- 20
pinter <- runif(ngenes)
lsobs <- -mean(log(pinter[genesetSize]))
system.time( lsperm <- mclapply(1:nresamp, function(i) -mean(log(sample(pinter, genesetSize)))) )
# system.time( lsperm <- mclapply(1:nresamp, function(i) -mean(log(sample(pinter, genesetSize))), mc.cores = 4) )
lsperm <- unlist(lsperm)
mean(lsperm > lsobs)
</syntaxhighlight>
* tempdir() is per session only. So the following will give the same directory.
:<syntaxhighlight lang='rsplus'>
parallel::mclapply(1:3, function(x) tempdir())
</syntaxhighlight>
* An example of using [https://statcompute.wordpress.com/2019/02/23/gradient-free-optimization-for-glmnet-parameters/ parallel::mcMap()].

Note
# Windows OS can not use mclapply(). The mclapply() implementation relies on forking and Windows does not support forking. mclapply from the parallel package is implemented as a serial function on Windows systems. The ''parallelsugar'' package was created based on the above idea.
# Another choice for Windows OS is to use parLapply() function in parallel package.
# [https://stackoverflow.com/questions/17196261/understanding-the-differences-between-mclapply-and-parlapply-in-r Understanding the differences between mclapply and parLapply in R] You don't have to worry about '''reproducing''' your environment on each of the cluster workers if mclapply() is used. <syntaxhighlight lang='rsplus'>
ncores <- as.integer( Sys.getenv('NUMBER_OF_PROCESSORS') )
cl <- makeCluster(getOption("cl.cores", ncores))
LLID2GOIDs2 <- parLapply(cl, rLLID, function(x) {
                                    library(org.Hs.eg.db); get("org.Hs.egGO")[[x]]} 
                        )
stopCluster(cl)
</syntaxhighlight>It does work. Cut the computing time from 100 sec to 29 sec on 4 cores.

== mclapply() vs parLapply() vs foreach() ==
* One of my examples shows mclapply() used much more memory than other methods
** ''mclapply''() uses 10.5GB while ''parLapply''() uses around 6.5GB and ''foreach''() (either doParallel or doSNOW backend) uses about 7GB (the numbers are obtained by monitoring the output from ''jobload'' command). NB: the memory used from the output of ''jobhist'' command in Biowulf is not reliable
** They are taking about the same amount of running time (19 minutes)
** mclapply() does not use ''clusterExport''() nor ''clusterEvalQ''() while others need to
** '''Bottomline: Do not use mclapply()''' if we want to save memory
** The code is available in my private bitbucket repository
* https://stackoverflow.com/questions/44806048/r-mclapply-vs-foreach
* [https://stackoverflow.com/a/11037275 mclapply does NOT copy the complete workspace contents to the child processes]
* [https://www.r-bloggers.com/trying-to-reduce-the-memory-overhead-when-using-mclapply/ Trying to reduce the memory overhead when using mclapply]. mclapply is the simplest (no need to export objects nor load packages) but it has a very high memory load.
* [https://www.uio.no/english/services/it/research/hpc/courses/r/2013-04_parallel_r.pdf Parallel programming in R] by Bjørn-Helge Mevik. parLapply() works on MPI for multiple machines.

== clusterApply() vs parLapply() and load-balancing versions ==
Look at the documentation [https://www.rdocumentation.org/packages/parallel/versions/3.6.2/topics/clusterApply clusterApply], it seems there is no difference in these two.

Note that there are load-balancing versions of these two: clusterApplyLB() and parLapplyLB().

== parallel vs doParallel package ==
doParallel is a '''foreach''' parallel adaptor for the '''parallel''' package.

doParallel provides a parallel backend for the '''%dopar%''' function using the parallel package.

doParallel depends on foreach, iterators, parallel and utils packages.

== parallelsugar package ==
* http://edustatistics.org/nathanvan/2015/10/14/parallelsugar-an-implementation-of-mclapply-for-windows/

If we load parallelsugar, the default implementation of parallel::mclapply, which used fork based clusters, will be overwritten by parallelsugar::mclapply, which is implemented with socket clusters. 

<syntaxhighlight lang='rsplus'>
library(parallel) 

system.time( mclapply(1:4, function(xx){ Sys.sleep(10) }) )
##    user  system elapsed 
##    0.00    0.00   40.06 

library(parallelsugar)
## 
## Attaching package: ‘parallelsugar’
## 
## The following object is masked from ‘package:parallel’:
## 
##     mclapply

system.time( mclapply(1:4, function(xx){ Sys.sleep(10) }) )
##    user  system elapsed 
##    0.04    0.08   12.98 
</syntaxhighlight>

== pbmcapply ==
[https://cran.r-project.org/web/packages/pbmcapply/index.html pbmcapply]: Tracking the Progress of Mc*pply with Progress Bar.

= snow package =
http://cran.r-project.org/web/packages/snow/index.html

Supported cluster types are "SOCK", "PVM", "MPI", and "NWS".

= [http://cran.r-project.org/web/packages/multicore/index.html multicore] package =
This package is removed from CRAN. 

Consider using package ‘parallel’ instead.

= [http://cran.r-project.org/web/packages/foreach/index.html foreach] package =

[https://blog.revolutionanalytics.com/2020/02/changes-in-the-foreach-package.html Changes in the foreach package] version 1.4.8, Feb 2020.

== Simple example ==
This package depends on one of the following 
* doParallel - Foreach parallel adaptor for the parallel package
* doSNOW - Foreach parallel adaptor for the snow package
* doMC - Foreach parallel adaptor for the multicore package. Used in [https://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html glmnet] vignette.
* doMPI - Foreach parallel adaptor for the Rmpi package
* doRedis - Foreach parallel adapter for the rredis package
as a backend.

<syntaxhighlight lang='rsplus'>
library(foreach)
library(doParallel)

m <- matrix(rnorm(9), 3, 3)

cl <- makeCluster(2, type = "SOCK")
registerDoParallel(cl) # register the parallel backend with the foreach package
foreach(i=1:nrow(m), .combine=rbind) %dopar%
  (m[i,] / mean(m[i,]))

stopCluster(cl)
</syntaxhighlight>

See also this post [http://blog.revolutionanalytics.com/2015/10/updates-to-the-foreach-package-and-its-friends.html Updates to the foreach package and its friends] on Oct 2015.

* [https://statcompute.wordpress.com/2015/12/13/calculate-leave-one-out-prediction-for-glm/ Cross validation in prediction for glm]
* [http://gforge.se/2015/02/how-to-go-parallel-in-r-basics-tips/#The_foreach_package How-to go parallel in R – basics + tips]

== combine list of lists ==
* .combine argument https://stackoverflow.com/questions/27279164/output-list-of-two-rbinded-data-frames-with-foreach-in-r
* [https://stackoverflow.com/questions/9519543/merge-two-lists-in-r Merge lists] by mapply() or base::Map()

<syntaxhighlight lang='rsplus'>
comb <- function(...) {
  mapply('cbind', ..., SIMPLIFY=FALSE)
}

library(foreach)
library(doParallel)

cl <- makeCluster(2)
registerDoParallel(cl) # register the parallel backend with the foreach package

m <- rbind(rep(1,3), rep(2,3))

# nrow(m) can represents number of permutations (2 in this toy example)
tmp <- foreach(i=1:nrow(m)) %dopar% {
  a <- m[i, ]
  b <- a * 10
  list(a, b)
}; tmp
# [[1]]
# [[1]][[1]]
# [1] 1 1 1
#
# [[1]][[2]]
# [1] 10 10 10
#
#
# [[2]]
# [[2]][[1]]
# [1] 2 2 2
#
# [[2]][[2]]
# [1] 20 20 20

foreach(i=1:nrow(m), .combine = "comb") %dopar% {
  a <- m[i,]
  b <- a * 10
  list(a, b)
}
# [[1]]
#      [,1] [,2]
# [1,]    1    2
# [2,]    1    2
# [3,]    1    2
#
# [[2]]
#      [,1] [,2]
# [1,]   10   20
# [2,]   10   20
# [3,]   10   20
stopCluster(cl)
</syntaxhighlight>

== Replacing double loops ==
* https://stackoverflow.com/questions/30927693/how-can-i-parallelize-a-double-for-loop-in-r
* http://www.exegetic.biz/blog/2013/08/the-wonders-of-foreach/
<syntaxhighlight lang='rsplus'>
library(foreach)
library(doParallel)

nc <- 4
nr <- 2

cores=detectCores()
cl <- makeCluster(cores[1]-1)
registerDoParallel(cl)
# set.seed(1234) # not work
# set.seed(1234, "L'Ecuyer-CMRG") # not work either
# library("doRNG")
# registerDoRNG(seed = 1985)    # not work with nested foreach
# Error in list(e1 = list(args = (1:nr)(), argnames = "i", evalenv = <environment>,  :
#  nested/conditional foreach loops are not supported yet.
m <- foreach (i = 1:nr, .combine='rbind') %:% # nesting operator
  foreach (j = 1:nc) %dopar% {
    rnorm(1, i*5, j) # code to parallelise
}
m
stopCluster(cl)
</syntaxhighlight>
Note that since the random seed (see the next session) does not work on nested loop, it is better to convert nested loop (two indices) to a single loop (one index).

== set seed and [https://cran.r-project.org/web/packages/doRNG/ doRNG] package ==
* [https://cran.r-project.org/web/packages/doRNG/vignettes/doRNG.pdf Vignette], [https://www.rdocumentation.org/packages/doRNG/versions/1.7.1/topics/doRNG-package Documentation]
* [http://michaeljkoontz.weebly.com/uploads/1/9/9/4/19940979/parallel.pdf#page=4 doRNG] package example
* [https://stackoverflow.com/questions/8358098/how-to-set-seed-for-random-simulations-with-foreach-and-domc-packages How to set seed for random simulations with foreach and doMC packages?]
* Use '''clusterSetRNGStream()''' from the parallel package; see [http://gforge.se/2015/02/how-to-go-parallel-in-r-basics-tips/ How-to go parallel in R – basics + tips]
* http://www.stat.colostate.edu/~scharfh/CSP_parallel/handouts/foreach_handout.html#random-numbers <syntaxhighlight lang='rsplus'>
library("doRNG") # doRNG does not need to be loaded after doParallel
library("doParallel")

cl <- makeCluster(2)
registerDoParallel(cl)

registerDoRNG(seed = 1234) # works for a single loop 
m1 <- foreach(i = 1:5, .combine = 'c') %dopar% rnorm(1)
registerDoRNG(seed = 1234)
m2 <- foreach(i = 1:5, .combine = 'c') %dopar% rnorm(1)
identical(m1, m2)
stopCluster(cl)

attr(m1, "rng") <- NULL # remove rng attribute
</syntaxhighlight>
* Another way to use the seed is to supply '''[https://www.rdocumentation.org/packages/doRNG/versions/1.7.1/topics/%25dorng%25 .options.RNG]''' in foreach() function. <syntaxhighlight lang='rsplus'>
r1 <- foreach(i=1:4, .options.RNG=1234) %dorng% { runif(1) }
</syntaxhighlight>

== Export libraries, variables, functions ==
* http://stat.ethz.ch/R-manual/R-devel/library/parallel/html/clusterApply.html
<syntaxhighlight lang='rsplus'>
clusterEvalQ(cl, {
  library(biospear)
  library(glmnet)
  library(survival)
})
clusterExport(cl, list("var1", "foo2"))
</syntaxhighlight>

If we use makeCluster(, "FORK"), then we don't need to do these.

== Summary the result ==
foreach returns the result in a list. For example, if each component is a matrix we can use

* Reduce("+", res)/length(res) # Reduce("+", res, na.rm = TRUE) not working
* apply(simplify2array(res), 1:2, mean, na.rm = TRUE)

to get the average of matrices over the list.

== Read a list files ==
[https://danielmarcelino.github.io/2018/reading-list-faster.html Reading List Faster With parallel, doParallel, and pbapply]

== Don't use too many cores: Error in unserialize(socklist<nowiki>[[n]]</nowiki>) ==
error reading from connection. See https://stackoverflow.com/a/24352032

== Why is foreach() sometimes slower than for? ==
See https://stackoverflow.com/a/10414280 or https://stackoverflow.com/q/16963808

* foreach is only advisable if you have relatively few rounds through very time consuming functions. 
* foreach() is best when the number of jobs does not hugely exceed the number of processors you will be using. 
<syntaxhighlight lang='rsplus'>
library(foreach)
library(parallel)
nCores <- future::availableCores()
cl <- makeCluster(nCores[[1]])
registerDoParallel(cl)

B<-10000

myFunc<-function() for(i in 1:B) sqrt(i)

myFunc2<-function() foreach(i = 1:B)  %do% sqrt(i)

myParFunc<-function() foreach(i = 1:B) %dopar% sqrt(i)

system.time(myFunc())
#   user  system elapsed 
#  0.001   0.002   0.004 
system.time(myFunc2())
#   user  system elapsed 
#  2.473   0.005   2.477 
system.time(myParFunc())
#   user  system elapsed 
#  3.560   0.242   3.829 

stopCluster(cl)
</syntaxhighlight>

Another example: reduce the probe sets corresponding to the same gene symbol by taking the probe set with the largest gene expression. It is not just the speed problem; it runs out of memory on a 44928x393 matrix with a PC has 8GB RAM. On a smaller data with 5000 genes and 100 arrays it took 0.071 seconds using split() + sapply() but it took 39 seconds using the foreach method.
<syntaxhighlight lang='rsplus'>
require(magrittr)
x <- matrix(rnorm(5000*100), nr=5000) # gene expression
y <- sample(3000, 5000, replace = TRUE) # gene symbol
rownames(x) <- names(y) <- 1:5000 # probe sets name

spt <- apply(x,1,mean,na.rm=TRUE) %>% split(., y)
ind <- sapply(spt, function(a) names(which.max(a)))
x.reduced <- x[ind, ]; rownames(x.reduced) <- names(spt)
</syntaxhighlight>

== Progress bar (doSNOW + foreach) ==
* [https://stackoverflow.com/a/40687156 How to show the progress of code in parallel computation in R?] It looks great.
{| class="wikitable"
! doParallel
! doSNOW
|-
| <pre>require(doParallel)
cores <- parallel::detectCores()
cl <- makeCluster(cores)
registerDoParallel(cl)
opts <- list(progress=progress)
result <- foreach(i=1:N) %dopar% {}
stopCluster(cl)</pre>
| <pre>require(doSNOW)
require(Kendall)
cores <- parallel::detectCores() 
cl <- makeSOCKcluster(cores) 
registerDoSNOW(cl) 
mydata <- matrix(rnorm(8000*500), ncol=500) 
pb <- txtProgressBar(min=1, max=8000, style=3)
progress <- function(n) setTxtProgressBar(pb, n)
opts <- list(progress=progress)
result <- foreach(i=1:8000, .packages="Kendall", 
          .options.snow=opts, .combine='rbind') %dopar% {
  abc <- MannKendall(mydata[i,])
  data.frame(tau=abc$tau, sl=abc$sl)
}
close(pb)
stopCluster(cl)</pre>
|}
* [https://blog.revolutionanalytics.com/2015/02/monitoring-progress-of-a-foreach-parallel-job.html Monitoring progress of a foreach parallel job] It also works on doSNOW but not doParallel backend.

= snowfall package =
http://www.imbi.uni-freiburg.de/parallel/docs/Reisensburg2009_TutParallelComputing_Knaus_Porzelius.pdf

= [http://cran.r-project.org/web/packages/Rmpi/index.html Rmpi] package =
Some examples/tutorials

* http://trac.nchc.org.tw/grid/wiki/R-MPI_Install
* http://www.arc.vt.edu/resources/software/r/index.php
* https://www.sharcnet.ca/help/index.php/Using_R_and_MPI
* http://math.acadiau.ca/ACMMaC/Rmpi/examples.html
* http://www.umbc.edu/hpcf/resources-tara/how-to-run-R.html
* [http://www.slideshare.net/bytemining/taking-r-to-the-limit-high-performance-computing-in-r-part-1-parallelization-la-r-users-group-727 Ryan Rosario]
* http://pj.freefaculty.org/guides/Rcourse/parallel-1/parallel-1.pdf
* * http://biowulf.nih.gov/apps/R.html

= OpenMP =
* [http://www.parallelr.com/r-and-openmp-boosting-compiled-code-on-multi-core-cpu-s/ R and openMP: boosting compiled code on multi-core cpu-s] from parallelr.com.

= [http://www.bioconductor.org/packages/release/bioc/html/BiocParallel.html BiocParallel] =
* [http://rpubs.com/seandavi/KallistoFromR Orchestrating a small, parallel, RNA-seq pre-processing workflow using R]

= [https://cran.r-project.org/web/packages/RcppParallel/index.html RcppParallel] =

= future & [https://cran.r-project.org/web/packages/future.apply/index.html future.apply] & doFuture packages =
* [https://alexioannides.com/2016/11/02/asynchronous-and-distributed-programming-in-r-with-the-future-package/ Asynchronous and Distributed Programming in R with the Future Package]
* [https://www.jottr.org/2018/06/23/future.apply_1.0.0/ Parallelize Any Base R Apply Function]
* [https://www.jottr.org/2019/01/11/parallelize-a-for-loop-by-rewriting-it-as-an-lapply-call/ Parallelize a For-Loop by Rewriting it as an Lapply Call]
* [https://blog.revolutionanalytics.com/2019/01/future-package.html Use foreach with HPC schedulers thanks to the future package]
* [https://www.jottr.org/2019/07/12/future-user2019-slides/ useR! 2019 Slides on Futures]
* [https://www.jottr.org/2019/11/09/resolved-launches-lazy-futures/ future 1.15.0 - Lazy Futures are Now Launched if Queried]
* [https://www.r-bloggers.com/rstudioconf-2020-slides-on-futures/ rstudio::conf 2020 Slides on Futures]

= Apache Spark =
* [http://files.meetup.com/3576292/Dubravko%20Dulic%20SparkR%20June%202016.pdf Introduction to Apache Spark]
* [https://therinspark.com/ Mastering Apache Spark with R] Javier Luraschi, Kevin Kuo, Edgar Ruiz (ebook)

= Microsoft R Server =
* [http://files.meetup.com/3576292/Stefan%20Cronjaeger%20R%20Server.pdf Microsoft R '''Server'''] (not Microsoft R Open)

= GPU =
* [http://www.parallelr.com/r-gpu-programming-for-all-with-gpur/ GPU Programming for All with ‘gpuR] from parallelr.com. The gpuR is available on [https://cran.r-project.org/web/packages/gpuR/index.html CRAN].
* [https://cran.r-project.org/web/packages/gputools/index.html gputools]

= Threads =
* [https://cran.r-project.org/web/packages/Rdsm/index.html Rdsm] package
* [https://random-remarks.net/2016/12/11/a-very-experimental-threading-in-r/ (A Very) Experimental Threading in R] and a post from [https://matloff.wordpress.com/2016/12/11/threading-in-r/ Mad Scientist]

= Benchmark =
* [http://rpsychologist.com/benchmark-parallel-sim Are parallel simulations in the cloud worth it? Benchmarking my MBP vs my Workstation vs Amazon EC2]

R functions to run timing
<syntaxhighlight lang='rsplus'>
# Method 1
system.time( invisible(rnorm(10000)))

# Method 2
btime <- Sys.time()
invisible(rnorm(10000))
Sys.time() - btime # difftime(Sys.time(), btime)
</syntaxhighlight>

* [https://www.r-bloggers.com/trying-to-reduce-the-memory-overhead-when-using-mclapply/ Trying to reduce the memory overhead when using mclapply]