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

\documentclass{beamer}

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

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

\begin{document}

<<setup, echo=FALSE>>=
options(width = 40)
@ 

\frame{\titlepage}

\begin{frame}{Motivation}
  Challenges
  \begin{itemize}
  \item Long calculations: bootstrap, MCMC, \ldots.
  \item Big data: genome-wide association studies, re-sequencing,
    \ldots.
  \item Long $\times$ big: \ldots
  \end{itemize}
  Solutions
  \begin{itemize}
  \item Avoid \R{} programming pitfalls -- \emph{very} significant
    benefits
  \item Large data management
  \item Parallel evaluation, especially `embarrassingly parallel' (not
    discussed in this course)
  \end{itemize}
\end{frame}

\section{Programming pitfalls}

\subsection{Pitfalls and solutions}

\begin{frame}[fragile]{Programming pitfalls: easy solutions}
  \begin{itemize}
  \item Input only required data
<<io-sketch, eval=FALSE, keep.source=TRUE>>=
colClasses <-
  c("NULL", "integer", "numeric", "NULL")
df <- read.table("myfile", colClasses=colClasses)
@ 
%% 
  \item Preallocate-and-fill, not copy-and-append
<<preallocate-and-fill-sketch, eval=FALSE, keep.source=TRUE>>=
result <- numeric(nrow(df))
for (i in seq_len(nrow(df)))
  result[[i]] <- some_calc(df[i,])
@ 

  \item Vectorized calculations, not iteration
<<iterative, eval=TRUE, keep.source=TRUE>>=
x <- runif(100000); x2 <- x^2
m <- matrix(x2, nrow=1000); y <- rowSums(m)
@ 

  \item Avoid unnecessary character creation operations, e.g.,
    \Rfunarg{USE.NAMES=FALSE} in \Rfunction{sapply}, \Rfunarg{use.names=FALSE}
    in \Rfunction{unlist}.
\end{itemize}
\end{frame}

\begin{frame}[fragile]{Programming pitfalls: moderate solutions}
  \begin{itemize}
  \item Use appropriate functions, often from specialized packages.
<<inappropriate-functions-2, eval=FALSE, keep.source=TRUE>>=
library(limma) # microarray linear models
fit <- lmFit(eSet, design)
@ 
  \item Identify appropriate algorithms, e.g., \Rfunction{\%in\%} is
    $O(N)$, whereas naive might be $O(N^2)$
<<algo-poly, keep.source=TRUE>>=
x <- 1:100; s <- sample(x, 10)
inS <- x %in% s
@ 
  \item Use C or Fortran code. Requires knowledge of other programming
    languages, and how to integrate these into \R{}
  \end{itemize}
\end{frame}

\subsection{Measuring performance}

\begin{frame}[fragile]{Measuring performance: timing}
  \begin{itemize}
  \item Use \Rfunction{system.time} to measure total evaluation time
    \begin{itemize}
    \item \Rfunarg{gcFirst=TRUE} for `garbage collection'
    \end{itemize}
  \item Use \Rfunction{replicate} to average over invocations
  \end{itemize}
<<system.time>>=
m <- matrix(runif(200000), 20000)
replicate(5, system.time(apply(m, 1, sum))[[1]])
replicate(5, system.time(rowSums(m))[[1]])
@ 
\begin{itemize}
\item Cautionary tale: \url{http://tinyurl.com/29bd6xv}
\end{itemize}
\end{frame}

\begin{frame}[fragile]{Measuring performance: comparison}
  \begin{itemize}
  \item \Rfunction{identical} and \Rfunction{all.equal} ensure that
    `optimizations' produce correct results!
  \end{itemize}
<<identical>>=
res1 <- apply(m, 1, sum)
res2 <- rowSums(m)
identical(res1, res2)
identical(c(1, -1), c(x=1, y=-1))
all.equal(c(1, -1), c(x=1, y=-1),
          check.attributes=FALSE)
@ 
\end{frame}

\begin{frame}[fragile]{Measuring execution time: \Rfunction{Rprof}}

<<Rprof, eval=FALSE, keep.source=TRUE>>=
tmpf = tempfile()
Rprof(tmpf)
res1 <- apply(m,  1, sum)
Rprof(NULL); summaryRprof(tmpf)
@ 
\begin{verbatim}
$by.self
         self.time self.pct total.time total.pct
"apply"       0.16       80       0.20       100
"FUN"         0.02       10       0.02        10
"lapply"      0.02       10       0.02        10
"unlist"      0.00        0       0.02        10

$by.total
         total.time total.pct self.time self.pct
"apply"        0.20       100      0.16       80
"FUN"          0.02        10      0.02       10
"lapply"       0.02        10      0.02       10
"unlist"       0.02        10      0.00        0
\end{verbatim}
\end{frame}

\begin{frame}[fragile]{Measuring memory use: \Rfunction{tracemem}}
  \begin{itemize}
  \item Enable memory profiling
  \end{itemize}
\begin{verbatim}
> ~/src/R-devel/configure --help
> ~/src/R-devel/configure --enable-memory-profiling
> make -j
\end{verbatim}
\begin{itemize}
\item Copy-on-change semantics
\end{itemize}
\begin{verbatim}
> x <- 1:10; tracemem(x)
[1] "<0x1b1a8f8>"
> y <- x       # no change, so no copy
> x[1] <- 2L   # x, y now differ, so copy
tracemem[0x1b1a8f8 -> 0x1b1a8a0]: 
\end{verbatim}
\end{frame}

\begin{frame}[fragile]{Measuring memory use: \Rfunction{tracemem}}
  \begin{itemize}
  \item Copying in \R{} functions
  \end{itemize}
\begin{verbatim}
> l <- list(a=1:10, b=1:10); tracemem(l$a)
[1] "<0x1131ce0>"
> df0 <- as.data.frame(l)
tracemem[0x1131ce0 -> 0x1131bd8]: eval as.data.frame.list as.data.frame 
tracemem[0x1131bd8 -> 0x1131a20]: data.frame eval eval as.data.frame.list as.data.frame 
tracemem[0x1131a20 -> 0x11318c0]: as.data.frame.integer as.data.frame data.frame eval eval as.data.frame.list as.data.frame 
> df1 <- data.frame(a=l$a, b=l$b)
tracemem[0x1131ce0 -> 0x11332c0]: data.frame 
tracemem[0x11332c0 -> 0x1133160]: as.data.frame.integer as.data.frame data.frame 
> identical(df0, df1)
[1] TRUE
\end{verbatim}
\end{frame}

\subsection{Case Study: GWAS}

\begin{frame}[fragile]{Case study: GWAS}
  \begin{itemize}
  \item Subset of genome-wide association study data
  \item Manage sample and SNP annotations in SQL 
  \item Manage SNP genotypes in netCDF
  \item Coordinate access using S4 classes
  \item Perform calculalations (sliding window composite linkage
    disequilibrium) in C.
  \item Organize as \R{} package \Rpackage{StudentGWAS}
  \end{itemize}
\end{frame}

\end{document}

%% Local Variables:
%% compile-command: "R CMD make -f devmakefile EfficientR.pdf"
%% End:
