---
title: Nesting `foreach` loops
author: Steve Weston
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{nested}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{utf8}
---

_Converted to RMarkdown by Hong Ooi_

## Introduction

```{r loadLibs, echo=FALSE, results="hide"}
library(foreach)
registerDoSEQ()
```

The `foreach` package provides a looping construct for executing
R code repeatedly.  It is similar to the standard `for` loop,
which makes it easy to convert a `for` loop to a `foreach`
loop.  Unlike many parallel programming packages for R, `foreach`
doesn't require the body of the `for` loop to be turned into a
function.  `foreach` differs from a `for` loop in that its
return is a list of values, whereas a `for` loop has no value and
uses side effects to convey its result.  Because of this,
`foreach` loops have a few advantages over `for` loops
when the purpose of the loop is to create a data structure such as a
vector, list, or matrix: First, there is less code duplication, and
hence, less chance for an error because the initialization of the vector
or matrix is unnecessary.  Second, a `foreach` loop may be easily
parallelized by changing only a single keyword.

## The nesting operator: `%:%`

An important feature of `foreach` is the `%:%` operator.
I call this the _nesting_ operator because it is used to create
nested `foreach` loops.  Like the `%do%` and
`%dopar%` operators, it is a binary operator, but it operates
on two `foreach` objects.  It also returns a `foreach`
object, which is essentially a special merger of its operands.

Let's say that we want to perform a Monte Carlo simulation using a
function called `sim`. (Remember that `sim` needs
to be rather compute intensive to be worth executing in parallel.)  The
`sim` function takes two arguments, and we want to call it with
all combinations of the values that are stored in the vectors
`avec` and `bvec`.  The following doubly-nested
`for` loop does that.  For testing purposes, the `sim`
function is defined to return $10 a + b$. (Of course, an
operation this trivial is not worth executing in parallel.)

```{r init1,echo=FALSE,results="hide"}
sim <- function(a, b) 10 * a + b
avec <- 1:2
bvec <- 1:4
```

```{r for1}
x <- matrix(0, length(avec), length(bvec))
for (j in 1:length(bvec)) {
  for (i in 1:length(avec)) {
    x[i,j] <- sim(avec[i], bvec[j])
  }
}
x
```

In this case, it makes sense to store the results in a matrix, so we
create one of the proper size called `x`, and assign the return
value of `sim` to the appropriate element of `x` each time
through the inner loop.

When using `foreach`, we don't create a matrix and assign values into
it.  Instead, the inner loop returns the columns of the result matrix as
vectors, which are combined in the outer loop into a matrix.
Here's how to do that using the `%:%` operator. Due to
operator precedence, you cannot put braces around the inner
`foreach` loop.

```{r foreach1}
x <-
  foreach(b=bvec, .combine='cbind') %:%
    foreach(a=avec, .combine='c') %do% {
      sim(a, b)
    }
x
```

This is structured very much like the nested `for` loop.
The outer `foreach` is iterating over the values in `bvec`,
passing them to the inner `foreach`, which iterates over the
values in `avec` for each value of `bvec`.  Thus, the `sim`
function is called in the same way in both cases.  The code is slightly
cleaner in this version, and has the advantage of being easily parallelized.

## Using `%:%` with `%dopar%`

When parallelizing nested `for` loops, there is always a question
of which loop to parallelize.  The standard advice is to parallelize the
outer loop.  This results in larger individual tasks, and larger tasks
can often be performed more efficiently than smaller tasks.  However, if
the outer loop doesn't have many iterations and the tasks are already
large, parallelizing the outer loop results in a small number of huge
tasks, which may not allow you to use all of your processors, and can
also result in load balancing problems.  You could parallelize an inner
loop instead, but that could be inefficient because you're repeatedly
waiting for all the results to be returned every time through the outer
loop.  And if the tasks and number of iterations vary in size, then it's
really hard to know which loop to parallelize.

But in our Monte Carlo example, all of the tasks are completely
independent of each other, and so they can all be executed in parallel.
You really want to think of the loops as specifying a single stream of
tasks.  You just need to be careful to process all of the results
correctly, depending on which iteration of the inner loop they came
from.

That is exactly what the `%:%` operator does: it turns multiple
`foreach` loops into a single loop.  That is why there is only
one `%do%` operator in the example above.  And when we
parallelize that nested `foreach` loop by changing the
`%do%` into a `%dopar%`, we are creating a single
stream of tasks that can all be executed in parallel:

```{r foreach2}
x <-
  foreach(b=bvec, .combine='cbind') %:%
    foreach(a=avec, .combine='c') %dopar% {
      sim(a, b)
    }
x
```

Of course, we'll actually only run as many tasks in parallel as we have
processors, but the parallel backend takes care of all that.  The point
is that the `%:%` operator makes it easy to specify the stream
of tasks to be executed, and the `.combine` argument to
`foreach` allows us to specify how the results should be processed.
The backend handles executing the tasks in parallel.

## Chunking tasks

Of course, there has to be a snag to this somewhere.  What if the tasks
are quite small, so that you really might want to execute the entire
inner loop as a single task?  Well, small tasks are a problem even for a
singly-nested loop.  The solution to this problem, whether you have a
single loop or nested loops, is to use _task chunking_.

Task chunking allows you to send multiple tasks to the workers at once.
This can be much more efficient, especially for short tasks.  Currently,
only the `doNWS` backend supports task
chunking.  Here's how it's done with `doNWS`:

```{r foreach3}
opts <- list(chunkSize=2)
x <-
  foreach(b=bvec, .combine='cbind', .options.nws=opts) %:%
    foreach(a=avec, .combine='c') %dopar% {
      sim(a, b)
    }
x
```

If you're not using `doNWS`, then this argument is ignored, which
allows you to write code that is backend-independent.  You can also
specify options for multiple backends, and only the option list that
matches the registered backend will be used.

It would be nice if the chunk size could be picked automatically, but I
haven't figured out a good, safe way to do that.  So for now, you need
to specify the chunk size manually.

The point is that by using the `%:%` operator, you can convert
a nested `for` loop to a nested `foreach` loop, use
`%dopar%` to run in parallel, and then tune the size of the
tasks using the `chunkSize` option so that they are big enough to be
executed efficiently, but not so big that they cause load balancing
problems.  You don't have to worry about which loop to parallelize,
because you're turning the nested loops into a single stream of tasks
that can all be executed in parallel by the parallel backend.

## Another example

Now let's imagine that the `sim` function returns a object that
includes an error estimate.  We want to return the result with the
lowest error for each value of b, along with the arguments that
generated that result.  Here's how that might be done with nested
`for` loops:

```{r init2, echo=FALSE, results="hide"}
sim <- function(a, b) {
  x <- 10 * a + b
  err <- abs(a - b)
  list(x=x, err=err)
}
```

```{r for2}
n <- length(bvec)
d <- data.frame(x=numeric(n), a=numeric(n), b=numeric(n), err=numeric(n))

for (j in 1:n) {
  err <- Inf
  best <- NULL
  for (i in 1:length(avec)) {
    obj <- sim(avec[i], bvec[j])
    if (obj$err < err) {
      err <- obj$err
      best <- data.frame(x=obj$x, a=avec[i], b=bvec[j], err=obj$err)
    }
  }
  d[j,] <- best
}
d
```

This is also quite simple to convert to `foreach`.  We just need
to supply the appropriate `.combine` functions.  For the outer
`foreach`, we can use the standard `rbind` function which can
be used with data frames.  For the inner `foreach`, we write a
function that compares two data frames, each with a single row,
returning the one with a smaller error estimate:

```{r innercombine}
comb <- function(d1, d2) if (d1$err < d2$err) d1 else d2
```

Now we specify it with the `.combine` argument to the inner
`foreach`:

```{r foreach4}
opts <- list(chunkSize=2)
d <-
  foreach(b=bvec, .combine='rbind', .options.nws=opts) %:%
    foreach(a=avec, .combine='comb', .inorder=FALSE) %dopar% {
      obj <- sim(a, b)
      data.frame(x=obj$x, a=a, b=b, err=obj$err)
    }
d
```

Note that since the order of the arguments to the `comb` function is
unimportant, I have set the `.inorder` argument to `FALSE`.
This reduces the number of results that need to be saved on the master
before they can be combined in case they are returned out of order.
But even with niceties such as parallelization, backend-specific
options, and the `.inorder` argument, the nested `foreach`
version is quite readable.

But what if we would like to return the indices into `avec` and
`bvec`, rather than the data itself?  A simple way to do that is to
create a couple of counting iterators that we pass to the
`foreach` functions:

```{r foreach5}
library(iterators)
opts <- list(chunkSize=2)
d <-
  foreach(b=bvec, j=icount(), .combine='rbind', .options.nws=opts) %:%
    foreach(a=avec, i=icount(), .combine='comb', .inorder=FALSE) %dopar% {
      obj <- sim(a, b)
      data.frame(x=obj$x, i=i, j=j, err=obj$err)
    }
d
```

Note that it's very important that the call
to icount is passed as the argument to `foreach`.  If the
iterators were created and passed to `foreach` using a variable,
for example, we would not get the desired effect.  This is not a bug or
a limitation, but an important aspect of the design of the
`foreach` function.

These new iterators are infinite iterators, but that's no problem since
we have `bvec` and `avec` to control the number of iterations of
the loops.  Making them infinite means we don't have to keep them in
sync with `bvec` and `avec`.

## Conclusion

Nested `for` loops are a common construct, and are often the most
time consuming part of R scripts, so they are prime candidates for
parallelization.  The usual approach is to parallelize the outer loop,
but as we've seen, that can lead to suboptimal performance due to an
imbalance between the size and the number of tasks.  By using
the `%:%` operator with `foreach`, and by using chunking
techniques, many of these problems can be overcome.  The resulting code
is often clearer and more readable than the original R code, since
`foreach` was designed to deal with exactly this kind of problem.