\documentclass{article}

\usepackage{BioconductorVignette}

\begin{document}

\SweaveOpts{eps=FALSE,width=4,height=4.5,include=FALSE, keep.source=TRUE}

\title{Efficient \R{} Programming}
\author{Martin Morgan \\
  Fred Hutchinson Cancer Research Center \\
  Seattle, WA, USA}
\date{July 30 2010}
\maketitle

\section{Pitfalls}

These brief exercises are meant to illustrate some common obstacles to
efficient \R{} programming. The idea is that you'll follow along with the
text, evaluating the instructions in your own \R{} session.

The basic scenario is a genome-wide association study. There are 1000
individuals. Case versus control status, gender, and age were recorded
for each, along with genotype at 10000 SNPs. The data is entirely
synthetic.

A suggestion from the lectures is to investigate packages that might
already efficiently implement the calculations you are trying to
perform. Check out the \Biocpkg{snpMatrix} package after doing these
exercises. How much of the following does this package make
unnecessary? Hint: almost all of it!

\subsection{Getting going}

To begin:
\begin{enumerate}
\item Start an \R{} session
\item Load the package for this conference.
<<load-pkg>>=
library(EfficientR)
@ 
\item Read several pre-defined functions in to the \R{} session
<<scripts>>=
fl <- system.file("script", "pitfalls.R", 
                  package="EfficientR")
source(fl)
@ 
\end{enumerate}

\subsection{Basic performance measurement and data I/O}

The first activities are meant to illustrate the use of
\Rfunction{system.time} to evaluate performance, \Rfunction{object.size} to
investigate how much memory an object uses, and \Rfunction{identical} and
\Rfunction{all.equal} to compare objects. We also look at how to more
effectively read data in to \R{}.

First, let's take a look at \Rfunction{fname0}, a file path to the GWAS
data set, and the function \Rfunction{f0}. The code below defines
\Rfunction{fname}; \Rfunction{f0} is defined in the \texttt{pitfalls.R}
file. Here are the definitions, and the result assigned to the
variable \Rcode{gwas1}
<<fname>>=
fname0 <- system.file("extdata", "gwas_2.csv", 
                      package="EfficientR")
f0
gwas1 <- f0(fname0)
@ 
%% 
This function simply reads in the file. Let's use \Rfunction{system.time}
to evaluate how long this takes
<<system.time-f0>>=
system.time(gwas0 <-  f0(fname0))
dim(gwas0)
@ 
%% 
Note the way in which the \Rfunction{<-} assignment to variable
\Rcode{gwas0} occurs in \Rfunction{system.time}. See the help page
\Rcode{?system.time} to understand the output; we're usually interested
in the \Rfunction{user.time}.

The function \Rfunction{f0} does its job, but perhaps we can improve on
it. For instance, suppose we were interested in a preliminary
investigation, and in particular reading just the case/control status,
sex, and age of the samples. The following shows two different
functions that allow us to read in just the information we're
interested in
<<f12>>=
f1
f2
@ 
%% 
Compare these functions with each other and \Rfunction{f0}, and with the
relevant help pages \Rcode{?read.csv}, \Rfunction{scan} to identify the key
steps for efficient data input. Now let's see how they perform in terms
of evaluation time:
<<system.time-f12>>=
system.time(gwas1 <- f1(fname0))
system.time(gwas2 <- f2(fname0))
@ 
%% 
Since we've read in just the data we're interested in, we should
have saved quite a bit of space
<<object.size>>=
object.size(gwas0)
object.size(gwas1)
@ 
%% 
The functions \Rfunction{f1} and \Rfunction{f2} are meant to return the same
result, just implemented in slightly different ways. We can verify
this with \Rfunction{identical}
<<identical-f12>>=
identical(gwas1, gwas2)    
@ 
%% 
Having read in the data, we might be interested in summarizing the
phenotypes we have, e.g.,
<<xtabs>>=
xtabs(~CaseControl + Sex, gwas2)
@ 

\subsection{Character manipulation}

Here we'll look at one aspect of \R{} that can sometimes have a
surprising performance penalty and sometimes tricky semantics:
character manipulation. We'll learn some additional techniques for
monitoring performance, as well as gain an appreciation for the
benefits of appropriate function choice.

Let's look at the genotypic data only, by dropping the first three
columns from the data frame; we'll take a peak at the first six rows
(\Rfunction{head}) of the first five columns of the genotype
information. (Pay no attention to the \Rfunction{names<-} call; it makes
the column names more distinct and is useful later).
<<gtype>>=
gtype <- gwas0[,-(1:3)]
names(gtype) <- paste("snp", names(gtype), sep="_")
head(gtype[,1:4])
@ 
%% 
Note that the rows have names (e.g., \Rcode{id\_1}).

A function \Rfunction{shuffle0} might come in handy if one were wanting to
randomize genotypes, and to return the randomized genotypes as a
matrix.
<<shuffle0>>=
shuffle0
@ 
%% 
We use \Rfunction{seed} to make the results reproducible across
invocations. \Rcode{sample(genotypes)} permutes the columns of the
\Rcode{gtype} data frame.  The \Rfunction{unlist} and \Rfunction{matrix} commands
are meant as a first attempt at creating a matrix from our permuted
data frame -- we `collapse' the sampled genotypes into a vector, and
then shape the vector into a matrix. Let's measure how long this
takes:
<<system.time-shuffl0>>=
system.time(s0 <- shuffle0(gtype))
@ 
%% 
Wow, that seems like a fairly long time for a relatively simple set of
operations. I wonder what's going on?
<<profile, keep.source=TRUE>>=
profFile <- tempfile()
Rprof(profFile)      # start gathering profile information
s0 <- shuffle0(gtype)
Rprof()              # stop
head(summaryRprof(profFile)$by.self)
@ 
%% 
A large fraction of the time is spent on the \Rfunction{unlist} function. A
little experimentation suggests what the problem might be:
<<unlist>>=
gsubset <- gtype[1:2, 1:3]
unlist(gsubset)
@ 
%% 
Notice that the value of \Rfunction{unlist} is a vector with names, and
that the names have been constructed to be unique. We can reference
the help page \Rcode{?unlist}, and arrive at a better solution
\Rfunction{shuffle1} that avoids creating names.
<<shuffle1>>=
shuffle1
@ 
%% 
This almost trivial change has a big influence on performance, without
changing the result:
<<system.time-shuffle1>>=
system.time(s1 <- shuffle1(gtype))
identical(s0, s1)
@ 
%% 
Finally, in \R{} it is common to be able to `cast' from one data
structure to another. \Rfunction{shuffle2} does this
<<shuffle2>>=
shuffle2
system.time(s2 <- shuffle2(gtype))
@ 
%% 
Note that the performance of \Rfunction{as.matrix} is comparable to our
\Rfunction{shuffle1}. Are the results the same?
<<identical-s12>>=
identical(s1, s2)
@ 
%% 
Oh oh! This doesn't look good. But maybe it's just that our results
\Rcode{s1} do not have \Rcode{dimnames}, whereas \Rcode{s2} might?
<<all.equal>>=
all.equal(s1, s2)
all.equal(s1, s2, check.attributes=FALSE)
@ 
%% 
(The result of the first call to \Rfunction{all.equal} is fairly
cryptic; this is, unfortunately, typical). The built-in function
\Rfunction{as.matrix} is performing as well as our version, and is
doing a better job of tracking important information (row and column
names) through the analysis.

\subsection{Preparing for an analysis}

To close this first exercise, let's set up a fairly typical but
statistically advanced GWAS-style analysis. I gloss over the details,
but we are doing a regression, where we try to explain the case versus
control classification in terms of sex, age, and a single SNP. We'd
ultimately be interested in SNPs that were most effective at
classifying a sample into case or control. Here's a function that
might perform the analysis of one SNP:
<<snp0>>=
snp0
snp0(1, gwas0)
@ 
%% 
These computations are expensive, especially when we remember that we
have (in real examples) hundreds of thousands of SNPs
<<system.time-snp>>=
system.time(result <- lapply(1:10, snp0, gwas0))
@ 
%% 
We'll see some approaches to getting better throughput in the next
section.

\section{Data I/O: Streaming}

We continue with the scenario of a genome-wide association
study. There are 1000 individuals. Case versus control status, gender,
and age were recorded for each, along with genotype at 10000 SNPs. The
data is entirely synthetic.

The following defines a function \Rfunction{fapply.csv} that (tries to!)
stream data from a file through a function, much like \Rfunction{lapply}
steams elements of a list through a function.
<<fapply>>=
readScript("fapply.R")
fl <- system.file("script", "fapply.R",
                  package="EfficientR")
source(fl)
@ 
%% 
By default, the function takes a file name, a function to be applied
to each chunk, and parameters influencing how the chunks are
processed. Let's give it a whirl, by reading the GWAS file in chunks
of 200 rows at a time. For each chunk, we'll calculate the fraction of
heterozygous loci. Here's our function
<<hetero>>=
hetero <- function(chunk, ...) {
  cat("starting chunk\n")
  rowSums(chunk[,-(1:3)]=="AB") / (ncol(chunk) - 3)
}
@ 
%% 
A common operation, both in stream processing and in distribution of
tasks for parallel evaluation, is to `reduce' the results from
different chunks / tasks into a single meaningful object. The
\Rfunarg{.reduce} argument to \Rfunction{fapply.csv} is meant to
be a function that performs the reduction. It expects a list, with
each element of the list the result of the function operating on a
chunk, e.g., each element being the result of \Rfunction{hetero}. As a
simple reduce function, we'll use \Rfunction{unlist}. The result
should be a vector of heterozygosities that we can manipulate or
visualize in the usual way. Let's give this a go\ldots
<<fapply-try>>=
res <- fapply.csv(fname0, hetero, row.names=1, .reduce=unlist)
length(res)
@ 
<<fapply-densityplot,fig=TRUE>>=
library(lattice)
print(densityplot(res, plot.points=FALSE, xlab="Heterozygosity", 
                  main="Stream"))
@ 
%% 
\begin{figure}
  \centering
  \includegraphics[height=!,width=.45\textwidth]{EfficientRProgrammingExercises-fapply-densityplot}
  \caption{Individual heterozygosity, `stream' processing}
\end{figure}
As exercises:
\begin{enumerate}
\item Explore different functions that might usefully be used in a
  streaming context. For instance, suppose the goal is to calculate
  heterozygosities for each SNP, rather than individual. What might
  the \Rfunarg{FUN} and \Rfunarg{.reduce} arguments to
  \Rfunction{fapply.csv} look like?
\item As an advanced exercise, download and investigate the
  \Rpkg{biglm} package for fitting linear models in a streaming
  fashion. When (what types of model and data) would the
  \Rpkg{biglm} package be particularly useful?
\item As another advanced exercise, consider how \Rfunction{fapply.csv}
  might be modified to work in parallel.
\end{enumerate}

\section{SQL}

Let's open the GWAS phenotype data from its SQLite data base\ldots
<<sqllite>>=
library(RSQLite)
db0 <- sub("csv$", "sql", fname0) # $
drv <- dbDriver("SQLite")
conn <- dbConnect(drv, dbname=db0)
dbListTables(conn)
@ 
%% 
\ldots and read all data into a data frame
<<sqlite-df>>=
q <- dbSendQuery(conn, "SELECT * FROM gwasPhenotypes")
df <- fetch(q, n=-1)
clear <- dbClearResult(q)
head(df)
@ 
%% 
\begin{enumerate}
\item How does the query result \Rcode{df} differ from the original
  data used to create the data base? Look carefully at the structure
  of the data frame, data values, and the class of columns in the data
  frame.  Does \Rfunction{all.equal} provide additional insight into
  differences or similarities? The original data can be retrieved with
<<orig-data>>=
gwasPhenotype <- f2(fname0)
@ 
\item What's the role of each of the commands starting with \Rcode{db}
  in the code above? use the help pages for assistance.
\item Two other ways of processing an entire table are
<<sql-table>>=
df1 <- dbGetQuery(conn, "SELECT * FROM gwasPhenotypes")
df2 <- dbReadTable(conn, "gwasPhenotypes")
@ 
\item SQL allows for more complicated queries. For instance evaluate
  the following and speculate on the meaning of the SQL components.
<<sql-which>>=
df <- dbGetQuery(conn, 
                 "SELECT age FROM gwasPhenotypes 
                  WHERE sex = 'F'")
@ 
%% 
\item SQL allows selection of several fields by providing a
  comma-delimiting string after 'select'. Can you formulate a SQL
  query to select the row names and \Rcode{Age} of all females?
\item SQL allows for functions to replace variable names as the
  'select' component. Verify (e.g., by comparison with the equivalent
  R statements on the entire data set) that the following returns the
  number of females over 40 in the data set.
<<sql-count>>=
dbGetQuery(conn,
           "SELECT COUNT(*) FROM gwasPhenotypes 
            WHERE age > 40 AND sex = 'F'")
@ 
\end{enumerate}
Don't forget to tidy up when done!
<<sql-tidy>>=
ok <- dbDisconnect(conn)
@ 

\section{NetCDF}

Turn now to the NetCDF file.
<<ncdf>>=
ncdf0 <- sub("csv$", "nc", fname0) #$
library(ncdf)
nc <- open.ncdf(ncdf0)
nc
@ 
%%
The \Rcode{nc} object contains useful summary information. It is in a
list-like structure. For instance, we can discover the number of
dimensions and their lengths with commands like
<<ncdf-discovery>>=
names(nc)
names(nc[["dim"]])
names(nc[["dim"]][["SNP"]])
nSnps <- nc[["dim"]][["SNP"]][["len"]]
@ 
\begin{enumerate}
\item What information can you `discover' about the \Rcode{Gentoype}
  variable?
\item Get the first 100 individuals and their SNPs with
<<ncdf-first-get>>=
g <- get.var.ncdf(nc, "Genotype", start=c(1, 1), count=c(10, nSnps))
@ 
\item From the examples on \Rcode{?get.var.ncdf}, is there another way
  to discover the dimensions of variables, and an easy way to retrieve
  all values in a dimension?
\item How would you retrieve samples 501 to 600, SNPs 1000 to 2000?
\item Would it be possible to retrieve all SNPs from females? What if
  the information from the \Rpkg{RSQLite} data were available?
  What kind of work-around might you imagine?
\end{enumerate}
%% 
Remember to tidy up when done!
<<ncdf-tidy>>=
ok <- close(nc)    
@ 

\section{\R{} on Clusters: \Rpkg{Rmpi}}

\textbf{These exercises require access to a cluster and so cannot be
  completed in this session; they may provide a useful resource.}
Commands necessary to run \R{} on a cluster tend to be very particular
to each system; the following uses \texttt{slurm} to manage cluster
jobs. A somewhat easier approach is to use multiple cores on a single
machine. The \Rpkg{multicore} and \Rpkg{foreach} packages
provide one example of this.

The basic scenario for these exercises is a genome-wide association
study. There are 1000 individuals. Case versus control status, gender,
and age were recorded for each, along with genotype at 10000 SNPs. The
data is entirely synthetic.

We've seen how to read in the data, and we're at the stage where we'd
like to fit generalized linear models to each SNP. We'll get to the
point where we can evaluate the model in parallel, and leave for an
`advanced exercise' steps necessary to make this actually a
computationally feasible endeavor.

\subsection{Getting going}

We start by becoming familiar with submitting batch jobs, and
exploring a common model for parallel evaluation in R.

Let's start by reviewing and submitting a batch job. Here's the
content of the file:
<<spawn>>=
readScript("spawn.R")
@ 
%% 
\begin{itemize}
\item Line 1 is a reminder about how we'll run this job. 

\item Line 2 loads the \Rpkg{Rmpi} library. \Rcode{mpi} (message passing
  interface) is a standard way to write parallel programs; \Rpkg{Rmpi}
  is a \emph{wrapper} around the \Rcode{mpi} libraries.

\item Line 3 consults the environment we are running in to find out
  how many computer nodes there are available for this process. We'll
  use one of the nodes for the `manager', the others for the `workers'
  (\Rpkg{Rmpi} uses the terminology `master' and `slave').

\item Line 4 starts the workers. Each worker starts on a particular
  machine, determined when we submit the batch job. The workers are
  assigned a \Rcode{rank} (0 for the manager, 1, 2, \ldots for the
  workers); we'll mention \Rcode{comm} latter. The worker starts a
  regular session of R, consuming just as many resources and with
  access to the same facilities as the manager.

\item Line 5 is an example of a \Rpkg{Rmpi} command. It evaluates the
  expression \Rfunction{mpi.comm.rank()} on each of the workers, and
  returns the result to the R session. The expression
  \Rfunction{mpi.comm.rank()} returns the rank of the machine on which it
  is evaluated; we should get one value for each worker.

\item Lines 6-7 define a function that takes a single argument, and
  returns that argument and the rank of the worker on which the
  function was evaluated.

\item Line 8 introduces a very useful way of using \Rpkg{Rmpi}:
  \Rfunction{mpi.parSapply} is like \Rfunction{sapply}, except that the
  calculation is distributed (as evenly as possible, by default)
  across the workers. So this line distributes \Rcode{1:(2*nWorkers)}
  `tasks' (i.e., determining \Rfunction{mpiRank}) across each worker. Other
  useful high-level functions include \Rfunction{mpi.parLapply} and
  \Rfunction{mpi.parReplicate}.

\item Lines 9 and 10 are `nice' ways to clean up, but are not actually
  required in the environment we're running in.
\end{itemize}
OK, this is how the job submission might look:
\begin{verbatim}
% salloc -N 4 mpirun -n 1 R --quiet -f spawn.R
\end{verbatim}
There are three parts to the job submission. First we use
\texttt{slurm}, the `job controller', to allocate 4 computer nodes for
our job. Then we use \texttt{mpirun} to launch, within the
\texttt{slurm} allocation, a single process. The process is \R{}, and
we launch \R{} in such a way that it is quiet (e.g., no start-up
message) and takes its input from the file \texttt{spawn.R}.
\texttt{slurm} responds by granting our allocation, launching
\texttt{mpirun} and then R, and then echoing the output back to the
screen:
\begin{verbatim}
salloc: Granted job allocation 239858
> # salloc -N 4 orterun -n 1 R -f spawn.R
> library(Rmpi)
> nWorkers <- mpi.universe.size() - 1L
> mpi.spawn.Rslaves()
        4 slaves are spawned successfully. 0 failed.
master (rank 0, comm 1) of size 5 is running on: hyraxB69
slave1 (rank 1, comm 1) of size 5 is running on: hyraxB69
slave2 (rank 2, comm 1) of size 5 is running on: hyraxB70
slave3 (rank 3, comm 1) of size 5 is running on: hyraxB71
slave4 (rank 4, comm 1) of size 5 is running on: hyraxB72
> mpi.remote.exec(mpi.comm.rank())
  X1 X2 X3 X4
1  1  2  3  4
> mpiRank <- function(i)
+   c(i=i, rank=mpi.comm.rank())
> mpi.parSapply(seq_len(2*nWorkers), mpiRank)
     [,1] [,2] [,3] [,4] [,5] [,6]
i       1    2    3    4    5    6
rank    1    2    2    3    3    4
> mpi.close.Rslaves()
[1] 1
> mpi.quit()
salloc: Relinquishing job allocation 239858
\end{verbatim}
We might capture the output by running the R script in `batch' mode
\begin{verbatim}
% salloc -N 4 mpirun -n 1 R CMD BATCH spawn.R
\end{verbatim}
which would generate a file \Rcode{spawn.Rout} in our current
directory. Another possibility is that we use standard Linux
redirection and piping to manipulate the output, e.g.,
\begin{verbatim}
% salloc -N 4 mpirun -n 1 R --vanilla < spawn.R | wc
\end{verbatim}

As exercises
\begin{enumerate}
\item Repeat the above actions in your session.
\item Try changing the number of nodes allocated by \texttt{salloc}
  (within reason!)
\item Investigate the \texttt{salloc} help page (\texttt{man salloc},
  using `f', `b' to page forward or backward, `q' to quit). Can you
  determine the difference between the -N and -n options?
\item Look at the scripts \texttt{random.R} and \texttt{simd.R} (use
  \Rcode{readScript("random.R")} and \Rcode{readScript("simd.R")}). Can
  you tell, from consulting with the class notes, what these scripts
  do?
\end{enumerate}

\subsection{Catching up}

In the first exercise, we set up a statistically advanced GWAS-style
analysis. We are doing a regression, where we try to explain the case
versus control classification in terms of sex, age, and a single SNP,
using a generalized linear model.  We can use this on a single node by
reading in the data and doing a number of analyses. Here's a simple
script:
<<snp0-script>>=
readScript("snp0.R")
@ 
%% 
We can invoke this on a single node in several ways; here are two
\begin{verbatim}
srun R -f snp0.R
salloc -N 1 mpirun -n 1 R -f snp0.R
\end{verbatim}
Now let's modify this script to support parallel evaluation. Here's a
first pass at the function:
<<snp2-script>>=
readScript("snp2.R")
@ 
%% 
The key features are:
\begin{itemize}
\item Lines 1 and 2: load \Rpkg{Rmpi} and spawn workers.
\item Lines 5-8: load the \Rcode{gwas} data set on each node, assigning
  to a global variable (with \Rfunction{<<-}).
\item Lines 11-15: define the function for fitting a single SNP,
  referencing the global variable. Return only the coefficients,
  rather than entire model.
\item Line 16: perform evaluation in parallel.
\end{itemize}
Note how we use the shared drives of the cluster to `distribute' the
data to each node, and restrict the amount of data that is returned to
the fitted coefficients (\Rfunction{glm} returns the original data along with
the fit, and so is large). One motivation for taking these steps is to
reduce the amount of explicit communication between nodes.

We can readily evaluate our script on 2, 5, or 10 nodes as, e.g.,
\begin{verbatim}
salloc -N 2 mpirun -n 1 R -f snp2.R
salloc -N 5 mpirun -n 1 R -f snp2.R
salloc -N 10 mpirun -n 1 R -f snp2.R
\end{verbatim}

\subsection{Looking forward}

The solution in \texttt{snp2.R} shows how to `get the job done', but
it is not very elegant, the code in \texttt{snp0.R} had to be modified
fairly extensively, and there are a number of redundant
calculations. As an advanced exercise,
\begin{enumerate}
\item The `ugliest' parts of the script are the use of a global
  variable \Rcode{gwas} to store the SNPs on each node, and the
  separation of the `data distribution' stage from the analysis
  stage. These are problems because they make the script more fragile
  and difficult to understand. Can you revise this script to address
  these issues? What about using a pattern like that in the
  \texttt{simd.R} (use \Rcode{readScript("simd.R")})?
\item For those with a strong statistical background. 

  One of the expensive parts of \Rfunction{glm} is the construction of a
  design matrix.  Note that most of the design matrix is constant
  across SNPs. Use this observation to construct a design matrix that
  excludes SNPs once, and explore how to modify this when a single SNP
  is added to the model.

  The \Rfunction{glm} algorithm depends on good starting values. Use this
  and the idea that the effect of individual SNPs is usually small to
  calculate starting values based on just the phenotypic data.

  Scanning all SNPs is likely a `first pass' at data analysis -- are
  there any SNPs that stand out as particularly interesting?
  \Rfunction{glm} is an iterative algorithm. What useful information
  can be obtained if the algorithm is restricted to just a single
  iteration?
\end{enumerate}

As an advanced exercise that integrates the data access and parallel
portions of the lab, here's a script that spawns two MPI nodes to
process (e.g., calculate heterozygosity per sample,
\Rfunction{nchetero}) different components of the data. Can you
identify the `reduce' functionality to transform the results of the
separate tasks into a single object? Can you transform the script
below into the `SIMD' style of MPI introduced last week, including
making the code scalable to different numbers of nodes?
<<ncdf-mpi, eval=FALSE>>=
nSamples <- 1000
nWorkers <- 2
## divide samples between workers; .splitIndices in Rmpi
library(Rmpi)
idx <- lapply(.splitIndices(nSamples, nWorkers), range)
## calculate individual heterozygosity for samples in idxElt
nchetero <- function(idxElt, ...)
{
    snpN <- nc[["dim"]][["SNP"]][["len"]]
    start <- c(idxElt[1], 1)
    count <- c(diff(idxElt) + 1L, snpN)
    d <- get.var.ncdf(nc, "Genotype", start=start, count=count)
    rowSums(d==2) / ncol(d)
}

mpi.spawn.Rslaves(nsl=nWorkers)
mpi.bcast.cmd(library(ncdf))
mpi.bcast.Robj2slave(ncdf0)
mpi.bcast.cmd(nc <- open.ncdf(ncdf0))
res <- mpi.parLapply(idx, nchetero)
mpi.bcast.cmd(close(nc))
ok <- mpi.close.Rslaves()
@ 
<<mpi-densityplot, eval=FALSE, fig=TRUE>>=
densityplot(unlist(res), plot.points=FALSE, xlab="Heterozygosity", 
            main="MPI and NetCDF")
@ 
%% 
Notice that this last example gains two speed-ups: I/O from use of
binary NetCDF files, and distributed computation across nodes.

\end{document}
