---
title: Using the `foreach` package
author: Steve Weston
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{foreach}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{utf8}
---

_Converted to RMarkdown by Hong Ooi_

## Introduction

One of R's most useful features is its interactive interpreter.  This
makes it very easy to learn and experiment with R.  It allows you to
use R like a calculator to perform arithmetic operations, display data
sets, generate plots, and create models.

Before too long, new R users will find a need to perform some
operation repeatedly.  Perhaps they want to run a simulation repeatedly
in order to find the distribution of the results.  Perhaps they need to
execute a function with a variety a different arguments passed to it.
Or maybe they need to create a model for many different data sets.

Repeated executions can be done manually, but it becomes quite
tedious to execute repeated operations, even with the use of command
line editing.  Fortunately, R is much more than an interactive
calculator.  It has its own built-in language that is intended to
automate tedious tasks, such as repeatedly executing R calculations.

R comes with various looping constructs that solve this problem.  The
`for` loop is one of the more common looping constructs, but
the `repeat` and `while` statements are also quite useful.
In addition, there is the family of "apply" functions, which includes 
`apply`, `lapply`, `sapply`, `eapply`,
`mapply`, `rapply`, and others. 

The `foreach` package provides a new looping construct for
executing R code repeatedly.  With the bewildering variety of existing
looping constructs, you may doubt that there is a need for yet another
construct.  The main reason for using the `foreach` package is
that it supports _parallel execution_, that is, it can execute those
repeated operations on multiple processors/cores on your computer, or on
multiple nodes of a cluster.  If each operation takes over a minute, and
you want to execute it hundreds of times, the overall runtime can take
hours.  But using `foreach`, that operation can be executed in
parallel on hundreds of processors on a cluster, reducing the execution
time back down to minutes.

But parallel execution is not the only reason for using the
`foreach` package.  There are other reasons that you might choose
to use it to execute quick executing operations, as we will see later in
the document.

## Getting Started

Let's take a look at a simple example use of the `foreach` package.
Assuming that you have the `foreach` package installed, you first
need to load it:

```{r loadLibs}
library(foreach)
```

Note that all of the packages that `foreach` depends on will be
loaded as well.

Now I can use `foreach` to execute the `sqrt` function
repeatedly, passing it the values 1 through 3, and returning the results
in a list, called `x`. (Of course, `sqrt` is a
vectorized function, so you would never really do this.  But later,
we'll see how to take advantage of vectorized functions with
`foreach`.)

```{r ex1}
x <- foreach(i=1:3) %do% sqrt(i)
x
```

This is a bit odd looking, because it looks vaguely like a `for`
loop, but is implemented using a binary operator, called
`%do%`.  Also, unlike a `for` loop, it returns a
value.  This is quite important.  The purpose of this statement is to
compute the list of results.  Generally, `foreach` with
`%do%` is used to execute an R expression repeatedly, and return
the results in some data structure or object, which is a list by
default.

You will note in the previous example that we used a variable `i` as
the argument to the `sqrt` function.  We specified the values of the 
`i` variable using a named argument to the `foreach` function.  We
could have called that variable anything we wanted, for example, `a`,
or `b`.  We could also specify other variables to be used in the R
expression, as in the following example:

```{r ex2}
x <- foreach(a=1:3, b=rep(10, 3)) %do% (a + b)
x
```

Note that parentheses are needed here.  We can also use braces:

```{r ex3}
x <- foreach(a=1:3, b=rep(10, 3)) %do% {
  a + b
}
x
```

We call `a` and `b` the _iteration variables_, since those are the
variables that are changing during the multiple executions.  Note that
we are iterating over them in parallel, that is, they are both changing
at the same time.  In this case, the same number of values are being
specified for both iteration variables, but that need not be the case.
If we only supplied two values for `b`, the result would be a list of
length two, even if we specified a thousand values for `a`:

```{r ex4}
x <- foreach(a=1:1000, b=rep(10, 2)) %do% {
  a + b
}
x
```

Note that you can put multiple statements between the braces, and you
can use assignment statements to save intermediate values of
computations.  However, if you use an assignment as a way of
communicating between the different executions of your loop, then your
code won't work correctly in parallel, which we will discuss later.

\section{The `.combine` Option}

So far, all of our examples have returned a list of results.  This is a
good default, since a list can contain any R object.  But sometimes
we'd like the results to be returned in a numeric vector, for example.
This can be done by using the `.combine` option to `foreach`:

```{r ex5}
x <- foreach(i=1:3, .combine='c') %do% exp(i)
x
```

The result is returned as a numeric vector, because the standard R `c`
function is being used to concatenate all the results.  Since the
`exp` function returns numeric values, concatenating them with
the `c` function will result in a numeric vector of length three.

What if the R expression returns a vector, and we want to combine those
vectors into a matrix?  One way to do that is with the `cbind` function:

```{r ex6}
x <- foreach(i=1:4, .combine='cbind') %do% rnorm(4)
x
```

This generates four vectors of four random numbers, and combines them by
column to produce a 4 by 4 matrix.

We can also use the `"+"` or `"*"` functions to combine our results:

```{r ex7}
x <- foreach(i=1:4, .combine='+') %do% rnorm(4)
x
```

You can also specify a user-written function to combine the results.
Here's an example that throws away the results:

```{r ex7.1}
cfun <- function(a, b) NULL
x <- foreach(i=1:4, .combine='cfun') %do% rnorm(4)
x
```

Note that this `cfun` function takes two arguments.  The
`foreach` function knows that the functions `c`, 
`cbind`, and `rbind` take many arguments, and 
will call them with up to 100 arguments (by default) in order to improve 
performance.  But if any
other function is specified (such as `"+"`), it assumes that it only
takes two arguments.  If the function does allow many arguments, you can
specify that using the `.multicombine` argument:

```{r ex7.2}
cfun <- function(...) NULL
x <- foreach(i=1:4, .combine='cfun', .multicombine=TRUE) %do% rnorm(4)
x
```

If you want the combine function to be called with no more than 10
arguments, you can specify that using the `.maxcombine` option:

```{r ex7.3}
cfun <- function(...) NULL
x <- foreach(i=1:4, .combine='cfun', .multicombine=TRUE, .maxcombine=10) %do% rnorm(4)
x
```

The `.inorder` option is used to specify whether the order in which the
arguments are combined is important.  The default value is
`TRUE`, but if the combine function is `"+"`, you could specify
`.inorder` to be `FALSE`.  Actually, this option is important 
only when executing the R expression in parallel, since results are always
computed in order when running sequentially.  This is not necessarily true when
executing in parallel, however.  In fact, if the expressions take very
different lengths of time to execute, the results could be returned in
any order.  Here's a contrived example, that executes the tasks in
parallel to demonstrate the difference.  The example uses the 
`Sys.sleep` function to cause the earlier tasks to take longer to execute:

```{r ex7.4}
foreach(i=4:1, .combine='c') %dopar% {
  Sys.sleep(3 * i)
  i
}
foreach(i=4:1, .combine='c', .inorder=FALSE) %dopar% {
  Sys.sleep(3 * i)
  i
}
```

The results of the first of these two examples is guaranteed to be the
vector `c(4, 3, 2, 1)`.  The second example will return the same values,
but they will probably be in a different order.

## Iterators

The values for the iteration variables don't have to be specified with
only vectors or lists.  They can be specified with an _iterator_, many
of which come with the `iterators` package.  An iterator is an
abstract source of data.  A vector isn't itself an iterator, but the
`foreach` function automatically creates an iterator from a
vector, list, matrix, or data frame, for example.  You can also create
an iterator from a file or a data base query, which are natural sources
of data.  The `iterators` package supplies a function called
`irnorm` which can return a specified number of random numbers
for each time it is called.  For example:

```{r ex8}
library(iterators)
x <- foreach(a=irnorm(4, count=4), .combine='cbind') %do% a
x
```

This becomes useful when dealing with large amounts of data.  Iterators
allow the data to be generated on-the-fly, as it is needed by your
operations, rather than requiring all of the data to be generated at the
beginning.

For example, let's say that we want to sum together a thousand random
vectors:

```{r ex9}
set.seed(123)
x <- foreach(a=irnorm(4, count=1000), .combine='+') %do% a
x
```

This uses very little memory, since it is equivalent to the following
`while` loop:

```{r ex10}
set.seed(123)
x <- numeric(4)
i <- 0
while (i < 1000) {
  x <- x + rnorm(4)
  i <- i + 1
}
x
```

This could have been done using the `icount` function, which
generates the values from one to 1000:

```{r ex11}
set.seed(123)
x <- foreach(icount(1000), .combine='+') %do% rnorm(4)
x
```

but sometimes it's preferable to generate the actual data with the
iterator (as we'll see later when we execute in parallel).

In addition to introducing the `icount` function from the
`iterators` package, the last example also used an unnamed
argument to the `foreach` function.  This can be useful when
we're not intending to generate variable values, but only controlling
the number of times that the R expression is executed.

There's a lot more that I could say about iterators, but for now,
let's move on to parallel execution.

## Parallel Execution

Although `foreach` can be a useful construct in its own right,
the real point of the `foreach` package is to do parallel computing.
To make any of the previous examples run in parallel, all you have to do
is to replace `%do%` with `%dopar%`.  But for the
kinds of quick running operations that we've been doing, there wouldn't
be much point to executing them in parallel.  Running many tiny tasks
in parallel will usually take more time to execute than running them
sequentially, and if it already runs fast, there's no motivation to make
it run faster anyway.  But if the operation that we're executing in
parallel takes a minute or longer, there starts to be some motivation.

### Parallel Random Forest

Let's take random forest as an example of an operation that can take
a while to execute.  Let's say our inputs are the matrix `x`, and the
factor `y`:

```{r ex12.data}
x <- matrix(runif(500), 100)
y <- gl(2, 50)
```

We've already loaded the `foreach` package, but we'll also need
to load the `randomForest` package:

```{r ex12.load}
library(randomForest)
```

If we want want to create a random forest model with a 1000 trees, and
our computer has four cores in it, we can split up the problem into four
pieces by executing the `randomForest` function four times, with
the `ntree` argument set to 250.  Of course, we have to combine
the resulting `randomForest` objects, but the
`randomForest` package comes with a function called
`combine` that does just that.

Let's do that, but first, we'll do the work sequentially:

```{r ex12.seq}
rf <- foreach(ntree=rep(250, 4), .combine=combine) %do%
  randomForest(x, y, ntree=ntree)
rf
```

To run this in parallel, we need to change `\%do\%`, but we also need to
use another `foreach` option called `.packages` to tell
the `foreach` package that the R expression needs to have the
`randomForest` package loaded in order to execute successfully.
Here's the parallel version:

```{r ex12.par}
rf <- foreach(ntree=rep(250, 4), .combine=combine, .packages='randomForest') %dopar%
  randomForest(x, y, ntree=ntree)
rf
```

If you've done any parallel computing, particularly on a cluster, you
may wonder why I didn't have to do anything special to handle `x` and
`y`.  The reason is that the `dopar` function noticed that
those variables were referenced, and that they were defined in the current
environment.  In that case `%dopar%` will automatically export
them to the parallel execution workers once, and use them for all of the
expression evaluations for that `foreach` execution.  That is
true for functions that are defined in the current environment as well,
but in this case, the function is defined in a package, so we had to
specify the package to load with the `.packages` option instead.

### Parallel Apply

Now let's take a look at how to make a parallel version of the standard
R `apply` function.  The `apply` function is written in R,
and although it's only about 100 lines of code, it's a bit difficult to
understand on a first reading.  However, it all really comes down two
`for` loops, the slightly more complicated of which looks like:

```{r ex13.orig}
applyKernel <- function(newX, FUN, d2, d.call, dn.call=NULL, ...) {
  ans <- vector("list", d2)
  for(i in 1:d2) {
    tmp <- FUN(array(newX[,i], d.call, dn.call), ...)
    if(!is.null(tmp)) ans[[i]] <- tmp
  }
  ans
}
applyKernel(matrix(1:16, 4), mean, 4, 4)
```

I've turned this into a function, because otherwise, R will complain
that I'm using `...` in an invalid context.

This could be executed using `foreach` as follows:

```{r ex13.first}
applyKernel <- function(newX, FUN, d2, d.call, dn.call=NULL, ...) {
  foreach(i=1:d2) %dopar%
    FUN(array(newX[,i], d.call, dn.call), ...)
}
applyKernel(matrix(1:16, 4), mean, 4, 4)
```

But this approach will cause the entire `newX` array to be sent
to each of the parallel execution workers.  Since each task needs only
one column of the array, we'd like to avoid this extra data
communication.

One way to solve this problem is to use an iterator that iterates over
the matrix by column:

```{r ex13.second}
applyKernel <- function(newX, FUN, d2, d.call, dn.call=NULL, ...) {
  foreach(x=iter(newX, by='col')) %dopar%
    FUN(array(x, d.call, dn.call), ...)
}
applyKernel(matrix(1:16, 4), mean, 4, 4)
```

Now we're only sending any given column of the matrix to one parallel
execution worker.  But it would be even more efficient if we sent the
matrix in bigger chunks.  To do that, we use a function called
`iblkcol` that returns an iterator that will return multiple columns
of the original matrix.  That means that the R expression will need to
execute the user's function once for every column in its submatrix.

```{r ex13.iter, results="hide"}
iblkcol <- function(a, chunks) {
  n <- ncol(a)
  i <- 1

  nextElem <- function() {
    if (chunks <= 0 || n <= 0) stop('StopIteration')
    m <- ceiling(n / chunks)
    r <- seq(i, length=m)
    i <<- i + m
    n <<- n - m
    chunks <<- chunks - 1
    a[,r, drop=FALSE]
  }

  structure(list(nextElem=nextElem), class=c('iblkcol', 'iter'))
}
nextElem.iblkcol <- function(obj) obj$nextElem()
```

```{r ex13.third}
applyKernel <- function(newX, FUN, d2, d.call, dn.call=NULL, ...) {
  foreach(x=iblkcol(newX, 3), .combine='c', .packages='foreach') %dopar% {
    foreach(i=1:ncol(x)) %do% FUN(array(x[,i], d.call, dn.call), ...)
  }
}
applyKernel(matrix(1:16, 4), mean, 4, 4)
```

Note the use of the `%do%` inside the `%dopar%` to
call the function on the columns of the submatrix `x`.  Now that
we're using `%do%` again, it makes sense for the iterator to be
an index into the matrix `x`, since `%do%` doesn't need to
copy `x` the way that `%dopar%` does.

## List Comprehensions

If you're familiar with the Python programming language, it may have
occurred to you that the `foreach` package provides something
that is not too different from Python's _list comprehensions_.
In fact, the `foreach` package also includes a function called
`when` which can prevent some of the evaluations from happening,
very much like the "if" clause in Python's list comprehensions.
For example, you could filter out negative values of an iterator using
`when` as follows:

```{r when}
x <- foreach(a=irnorm(1, count=10), .combine='c') %:% when(a >= 0) %do% sqrt(a)
x
```

I won't say much on this topic, but I can't help showing how
`foreach` with `when` can be used to write a simple quick
sort function, in the classic Haskell fashion:

```{r qsort}
qsort <- function(x) {
  n <- length(x)
  if (n == 0) {
    x
  } else {
    p <- sample(n, 1)
    smaller <- foreach(y=x[-p], .combine=c) %:% when(y <= x[p]) %do% y
    larger  <- foreach(y=x[-p], .combine=c) %:% when(y >  x[p]) %do% y
    c(qsort(smaller), x[p], qsort(larger))
  }
}

qsort(runif(12))
```

Not that I recommend this over the standard R `sort` function.
But it's a pretty interesting example use of `foreach`.

## Conclusion

Much of parallel computing comes to doing three things: splitting the
problem into pieces, executing the pieces in parallel, and combining the
results back together.  Using the `foreach` package, the
iterators help you to split the problem into pieces, the
`%dopar%` function executes the pieces in parallel, and the
specified `.combine` function puts the results back together.
We've demonstrated how simple things can be done in parallel quite
easily using the `foreach` package, and given some ideas about
how more complex problems can be solved.  But it's a fairly new package,
and we will continue to work on ways of making it a more powerful system
for doing parallel computing.