%\VignetteIndexEntry{Efficient R Programming -- labs}
%\VignettePackage{AdvancedR2011}

\documentclass{article}

\usepackage{mathptmx}

\usepackage{BioconductorSlides}
\SweaveOpts{keep.source=TRUE, eval=TRUE}

\title{Efficient \R{} Programming}
\author{Martin Morgan \\
 Fred Hutchinson Cancer Research Center}
\date{17-18 February, 2011}

\begin{document}

\SweaveOpts{eps=FALSE,width=4,height=4.5,include=FALSE, eval=TRUE}

\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 $\approx$100000 SNPs. The data is
entirely synthetic.

\subsection{Getting going}

To begin:
\begin{enumerate}
\item Start an \R{} session, and use \Rfunction{sessionInfo()} to
  confirm that you have version 2.12.1.
\item Install the \Rpackage{AdvancedR2011Data} package from the thumb
  drive distributed at the start of the class.
\item Install the course package with
\begin{verbatim}
  source("http://bioconductor.org/course-packages/courseInstall.R")
  courseInstall("AdvancedR2011")
\end{verbatim}
\item Ensure that the installation has been successful by loading the
  \Rpackage{AdvancedR2011} package
<<library>>=
library(AdvancedR2011)
@ 
\item Read several pre-defined functions in to the R session
<<scripts>>=
fl <- system.file("script", "efficient.R", 
                  package="AdvancedR2011")
source(fl)
@ 
\end{enumerate}

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

The first activities are meant to illustrate the use of
\Rcode{system.time} to evaluate performance, \Rcode{object.size} to
investigate how much memory an object uses, and \Rcode{identical} and
\Rcode{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 \Rcode{fname0}, a file path to the GWAS
data, and the function \Rcode{f0}. Both \Rcode{fname} and \Rcode{f0}
are defined in the \texttt{efficient.R} file. Here are the
definitions, and the result assigned to the variable \Rcode{gwas0}; we
read in just 3 rows of data, to keep the exercise managable. Let's use
\Rcode{system.time} to evaluate how long this takes
<<fname>>=
fname0
f0
system.time(gwas0 <- f0(fname0, nrows=3))
dim(gwas0)
@ 
%% 
Note the way in which the \Rcode{<-} assignment to variable
\Rcode{gwas0} occurs in \Rcode{system.time}. See the help page
\Rcode{?system.time} to understand the output; we're usually
interested in the \Rcode{user.time}.

The function \Rcode{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 first 100 SNPs. 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 \Rcode{f0}, and with the
relevant help pages \Rcode{?read.csv}, \Rcode{scan} to identify the
key steps for efficient data input. Now lets see how they perform in
terms of evaluation time:
<<system.time-f12>>=
system.time(gwas1 <- f1(fname0, nrows=3))
system.time(gwas2 <- f2(fname0, nrows=3))
@ 
%% 
Since we've read in the 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 \Rcode{f1} and \Rcode{f2} are meant to return the same
result, just implemented in slightly different ways. We can verify
this with \Rcode{identical}
<<identical-f12>>=
identical(gwas1, 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 the genotypic data only, by dropping the first three
columns from the data frame; we'll take a peak at the first six rows
(\Rcode{head}) of the first five columns of the genotype
information.
<<gtype>>=
gtype <- f2(fname0, keep=1:10000, nrows=-1)
head(gtype[,1:5])
@ 
%% 
Note that the columns have names (e.g., \Rcode{id\_V1}).

A function \Rcode{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 \Rcode{seed} to make the results reproducible across
invocations. \Rcode{sample(genotypes)} permutes the columns of the
\Rcode{gtype} data frame.  The \Rcode{unlist} and \Rcode{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(NULL)          # stop
head(summaryRprof(profFile)$by.self)
@ 
%% 
A large fraction of the time is spent on the \Rcode{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 \Rcode{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
\Rcode{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{} its common to be able to `cast' from one data
structure to another. \Rcode{shuffle2} does this
<<shuffle2>>=
shuffle2
system.time(s2 <- shuffle2(gtype))
@ 
%% 
Note that the performance of \Rcode{as.matrix} is comparable to our
\Rcode{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 \Rcode{all.equal} is fairly cryptic;
this is, unfortunately, typical.) The built-in function
\Rcode{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.

\section{Data I/O: Streaming}

It can be difficult to fit large data into memory, and it is
increasingly necessary to think of `streaming' algorithms that process
a small `chunk' of data at a time.  The following defines a function
\Rfunction{.fapply} that (tries to!) stream data from a file through a
function, much like \Rfunction{lapply} streams elements of a list
through a function.
<<fapply>>=
.fapply
@ 
%%
The function takes an open file connecton, a function to be applied to
each chunk, and parameters influencing how the chunks are input and
processed. Let's give it a whirl, by reading the GWAS file in chunks
of 100 rows at a time. For each chunk, we'll calculate the number of
individuals in the chunk, and the heterozygotes for each SNP. Here's
our function
<<fapply-FUN>>=
heterozygosity <- function(chunk, ...)
    list(N=nrow(chunk), Het=colSums(chunk == 2))
@ 
%%
The \Rfunarg{.get} argument to \Rfunction{.fapply} is meant to retrieve
a chunk of data, and takes as its single argument an open \R{}
connection. Here's our \Rfunction{get}:
<<fapply-get>>=
get <- function(con)
    f2(con, nrows=100, ncols=.ncols, keep=1:100)
@ 
%% 
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} 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{heterozygosity}. The
result should be a vector of heterozygosities that we can manipulate
or visualize in the usual way. Here's our simple reduce function:
<<fapply-reduce>>=
reduce <- function(lst)
{
    n <- sum(sapply(lst, "[[", "N"))
    het <- Reduce("+", lapply(lst, "[[", "Het")) 
    browser()
    het / n
}
@ 
Let's give this a go\ldots
<<fapply>>=
con <- file(fname0, open="r")
res <- .fapply(con, heterozygosity, .get=get, .reduce=reduce)
close(con)
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]{EfficientR-lab-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 individual, rather than each SNP. What
  might the \Rfunarg{FUN} and \Rfunarg{.reduce} arguments to
  \Rfunction{.fapply} look like?
\item As an advanced exercise, consider how \Rfunction{.fapply} might
  be modified to work in parallel, for instance using the
  \Rpackage{multicore} or \Rpackage{Rmpi} packages.
\item The \Rfunarg{.reduce} argument takes the complete list of
  results, and reduces it. This will not be a good strategy if the
  result of \Rfunarg{FUN} is itself large. Revise \Rfunction{.fapply}
  so that \Rfunarg{.reduce} is a function that takes as its first
  argument the current reduction, and as its second argument the
  current chunk with \Rfunarg{FUN} already applied.
\end{enumerate}

\end{document}
