
%\VignetteIndexEntry{matter: Rapid prototyping with data on disk}
%\VignetteKeyword{Infrastructure}

\documentclass[a4paper]{article}
\usepackage{caption}
\usepackage{subcaption}


<<eval=TRUE, echo=FALSE, results=tex>>=
BiocStyle::latex()
@

\title{\Rpackage{matter}: Rapid prototyping with data on disk}

\author{Kylie A. Bemis}

\begin{document}
\SweaveOpts{concordance=TRUE}

\maketitle

\tableofcontents

\section{Introduction}

\Rpackage{matter} is designed for rapid prototyping of new statistical methods when working with larger-than-memory datasets on disk. Unlike related packages \Rpackage{bigmemory} \cite{bigmemory} and \Rpackage{ff} \cite{ff}, which also work with file-backed larger-than-memory datasets, \Rpackage{matter} aims to offer strict control over memory and maximum flexibility with on-disk data structures, so it can be easily adapted to domain-specific file formats, including user-customized file formats.

\section{Installation}

\Rpackage{matter} can be installed from Bioconductor using the following commands in \R{}.

<<eval=FALSE>>=
source("https://bioconductor.org/biocLite.R")
biocLite("matter")
@

\section{Basic use and data manipulation}

<<echo=FALSE,results=hide>>=
options(width=72)
library(matter)
data(matter_ex)
@ 

\Rpackage{matter} matrices and vectors can be initialized similarly to ordinary \R{} matrices. When no file is given, a new temporary file is created in the default temporary file directory, which will be cleaned up later by either \R{} or the operating system.

Here, we initialize a \Rpackage{matter} matrix with 10 rows and 5 columns. The resulting object is a subclass of the \Robject{matter} class, and stores file metadata that gives the location of the data on disk. In many cases, it can be treated as an ordinary \R{} matrix.

<<>>=
x <- matter_mat(data=1:50, nrow=10, ncol=5)
x
x[]
@

As seen above, this is a small toy example in which the in-memory metadata actually takes up more space than the size of the data stored on disk. For much larger datasets, the in-memory metadata will be a small fraction of the total size of the dataset on disk.

\Rpackage{matter}'s matrices and vectors can be indexed into like ordinary \R{} matrices and vectors.

<<>>=
x[1:4,]
x[,3:4]
@

We can assign names to \Robject{matter\_vec} vectors and row and column names to \Robject{matter\_mat} matrices.

<<>>=
rownames(x) <- 1:10
colnames(x) <- letters[1:5]
x[]
@

\Rpackage{matter} provides methods for calculating summary statistics for its vectors and matrices, including some methods that do not exist in base \R{}, such as \Robject{colVars}, which uses a memory-efficient running variance calculation  that is accurate for large floating-point datasets \cite{Welford:1962tw}.

<<>>=
colSums(x)
colSums(x[])

colVars(x)
apply(x, 2, var)
@

One of the major advantages of the flexibility of \Rpackage{matter} is being able to treat data from multiple files as a single dataset. This is particularly useful if analysing data from a domain where each sample in an experiment generates large files, such as high-resolution, high-throughput mass spectrometry imaging.

Below, we create a second matrix, and show its data is stored in a separate file. We then combine the matrices, and the result can be treated as a single matrix, despite originating from multiple files. Combinging the matrices does not create new data or change the existing data on disk.

<<>>=
y <- matter_mat(data=51:100, nrow=10, ncol=5)

paths(x)
paths(y)

z <- cbind(x, y)
z
z[]
@

Note that matrices in \Rpackage{matter} are either stored in a column-major or a row-major format. The default is to use the column-major format, as \R{} does. Column-major matrices are optimized for fast column-access, and assume that each column is stored contiguously or mostly-contiguously on disk. Conversely, row-major matrices are optimized for fast row-access, and make the same assumption for rows.

Since \Rpackage{matter} does support both column-major and row-major formats, transposing a matrix is a trivial operation in \Rpackage{matter} that only needs to change the matrix metadata, and doesn't touch the data on disk.

<<>>=
t(x)

rbind(t(x), t(y))
@

Note that this is equivalent to \verb|t(cbind(x, y))|.

Below, we inspect the metadata associated with the different columns of \Robject{x} using the \verb|adata| method.

<<>>=
adata(x)
@

Note that each column has a byte offset and an extent (i.e., length) associated with it.

Now we show how to create a \Robject{matter\_mat} matrix for an pre-existing file. We will point the new matrix to the bottom half of \Robject{x}.

<<>>=
xsub <- matter_mat(offset=c(40, 120, 200, 280, 360),
            extent=rep(5,5), paths=paths(x))
x[6:10,]
xsub[]
@

It is possible to build \Robject{matter} objects from nearly any possible combination of files and locations within files. It is even possible to build up a matrix from vectors, which we do below.

<<>>=
x2 <- matter_vec(offset=80, extent=10, paths=paths(x))
y3 <- matter_vec(offset=160, extent=10, paths=paths(y))
cbind(x2, y3)[]
cbind(x[,2], y[,3])
@

This is a quick and easy way to build a dataset from many files where each column of the dataset is stored in a separate file. Even if the resulting matrix would fit into memory, using \Rpackage{matter} can be a tidy, efficient way of reading complex binary data from multiple files into \R{}.


\section{Linear regression for on-disk datasets}

\Rpackage{matter} is designed to provide a statistical computing environment for larger-than-memory datasets on disk. To facilitate this, \Rpackage{matter} provides a method for fitting of linear models for \Robject{matter} matrices through the \Rpackage{biglm} package \cite{biglm}. \Rpackage{matter} provides a wrapper for \Rpackage{biglm}'s \Robject{bigglm} function that works with \Robject{matter\_mat} matrices, which we demonstrate below.

First, we simulate some data appropriate for linear regression.

<<setup-bigglm>>=
set.seed(81216)
n <- 15e6
p <- 9

b <- runif(p)
names(b) <- paste0("x", 1:p)

data <- matter_mat(nrow=n, ncol=p + 1)

colnames(data) <- c(names(b), "y")

data[,p + 1] <- rnorm(n)
for ( i in 1:p ) {
  xi <- rnorm(n)
  data[,i] <- xi
  data[,p + 1] <- data[,p + 1] + xi * b[i]
}

data

head(data)
@

This creates a 1.2 GB dataset on disk, but less than 20 KB of metadata is stored in memory.

Now we calculate some statistical summaries using \Rpackage{matter}'s \Robject{apply} method for \Robject{matter\_mat} matrices.

<<summary-stats>>=
apply(data, 2, mean)

apply(data, 2, var)
@

We could also have used \Robject{colMeans} and \Robject{colVars}.

Now we fit the linear model to the data using the \Robject{bigglm} method for \Robject{matter\_mat} matrices. Note that it requires a formula, and (unfortunately) it does not allow \verb|y ~ .|, so all variables must be stated explicitly.

%%%%%------- offline analysis --------%%%%%%

% <<fit-bigglm-offline>>=
% matter:::profile({
%   fm <- as.formula(paste0("y ~ ", paste(names(b), collapse=" + ")))
%   bigglm.out <- bigglm(fm, data=data, chunksize=10000)
% })
%    # start (MB)   finish (MB) max used (MB) overhead (MB)    time (sec) 
%     # 206.400       207.200       702.600       495.400        46.264 
%     # 206.400       207.200       702.600       495.400        45.717
%     # 206.400       207.200       702.600       495.400        44.197 
%     # 206.400       207.200       702.600       495.400        45.988

% data.m <- data[]
% data.df <- as.data.frame(data.m)

% rm(data.m)
% gc()

% matter:::profile({
%   lm.out <- lm(fm, data=data.df)
% })
%    # start (MB)   finish (MB) max used (MB) overhead (MB)    time (sec) 
%    # 1406.500      5660.600      7100.900      1440.300        36.226 
%      # 1406.500      5660.600      7100.900      1440.300        33.151
%      # 1406.500      5660.600      7100.900      1440.300        31.838

% matter:::profile({
%   biglm.out <- bigglm(fm, data=data.df, chunksize=10000)
% })
%    # start (MB)   finish (MB) max used (MB) overhead (MB)    time (sec) 
%    # 1406.500      1410.700      2942.100      1531.400       154.417
%    # 1406.500      1410.700      2942.100      1531.400       157.163
%    # 1406.500      1410.700      2942.100      1531.400       158.186

% rm(data.df)

% require(biganalytics)
% data.m <- data[]
% data.bm <- as.big.matrix(data.m)

% rm(data.m)
% gc()

% matter:::profile({
%   biglm.bm.out <- biglm.big.matrix(fm, data=data.bm, chunksize=10000)
% })
%    # start (MB)   finish (MB) max used (MB) overhead (MB)    time (sec) 
%    # 211.100       211.600      1378.900      1167.300        20.848 
%    # 211.100       211.600      1378.900      1167.300        21.267
%    # 211.100       211.600      1378.900      1167.300        21.621 

% require(ffbase)    
% data.ff <- ff(vmode="double", dim=dim(data.m))
% data.ff[] <- data.m
% data.ff <- as.ffdf(data.ff)
% colnames(data.ff) <- colnames(data)

% rm(data)
% rm(data.m)
% gc()

% matter:::profile({
%   biglm.ff.out <- bigglm(fm, data=data.ff, chunksize=10000)
% })
%    # start (MB)   finish (MB) max used (MB) overhead (MB)    time (sec) 
%    #    209.700       218.700       647.400       428.700        57.182 
%        # 209.700       218.700       647.400       428.700        59.032 
%        # 209.70        218.70        647.40        428.70         58.52
%        # 209.700       218.700       647.400       428.700        58.561
%        # 209.70        218.70        647.40        428.70         58.41
% @

%%%%%------- online analysis --------%%%%%%

<<fit-bigglm, eval=FALSE>>=
fm <- as.formula(paste0("y ~ ", paste(names(b), collapse=" + ")))
bigglm.out <- bigglm(fm, data=data, chunksize=10000)
@

<<coef>>=
summary(bigglm.out)

cbind(coef(bigglm.out)[-1], b)
@

On a 2012 retina MacBook Pro with 2.6 GHz Intel CPU, 16 GB RAM, and 500 GB SSD, fitting the linear model takes 41 seconds and uses an additional 257 MB of memory overhead. The max amount of memory used while fitting the model was only 499 MB, for the 1.2 GB dataset. This memory usage can be controlled further by using the \verb|chunksize| argument in \verb|bigglm| or by specifying a different \verb|chunksize| for the \Robject{matter} object.


\section{Principal components analysis for on-disk datasets}

Because \Rpackage{matter} provides basic linear algebra for on-disk \Robject{matter\_mat} matrices with in-memory R matrices, it opens up the possibility for the use of many iterative statistical methods which can operate on only small portions of the data at a time.

For example, \Robject{matter\_mat} matrices are compatible with the \Rpackage{irlba} package, which performs efficient, bounded-memory singular value decomposition (SVD) of matrices, and which can therefore be used for efficient principal components analysis (PCA) of large datasets \cite{irlba}.

For convenience, \Rpackage{matter} provides a \verb|prcomp| method for performing PCA on \Robject{matter\_mat} matrices using the \verb|irlba| method from the \Rpackage{irlba} package, as demonstrated below.

First, we simulate some data appropriate for principal components analysis.

<<setup-prcomp>>=
set.seed(81216)
n <- 15e5
p <- 100

data <- matter_mat(nrow=n, ncol=p)

for ( i in 1:10 )
  data[,i] <- (1:n)/n + rnorm(n)
for ( i in 11:20 )
  data[,i] <- (n:1)/n + rnorm(n)
for ( i in 21:p )
  data[,i] <- rnorm(n)

data
@

This again creates a 1.2 GB dataset on disk, but only 160 KB of metadata is stored in memory. More metadata is stored compared to the previous example, because the matrix has more columns, and it is stored as a column-major \Robject{matter\_matc} matrix with independent metadata for each column.

Note that, in the simulated data, only the first twenty variables show systematic variation, with the first ten variables varying distinctly from the next ten variables.

First we calculate the variance for each column.

<<plot-var>>=
var.out <- colVars(data)
plot(var.out, type='h', ylab="Variance")
@

This takes only 7 seconds and uses less than 30 KB of additional memory. The maximum amount of memory used while calculating the variance for all columns of the 1.2 GB dataset is only 27 MB.

Note that the \Robject{irlba} function has an optional argument \verb|mult| which allows specification of a custom matrix multiplication method, for use with packages such as \Rpackage{bigmemory} and \Rpackage{ff}. This is especially useful since it allows a \verb|transpose=TRUE| argument, so that the identity \verb|t(t(B) %*% A)| can be used in place of \verb|t(A) %*% B)| when transposition is an expensive operation. However, this is not necessary for \Rpackage{matter}, since transposition is a trivial operation for \Robject{matter\_mat} matrices.

Implicit scaling and centering is supported in \Rpackage{matter}, so the data can be scaled and centered without changing the data on disk or the additional overhead of storing a scaled version of the data.

Unlike the default \verb|prcomp| method for ordinary \R{} matrices, this version supports an argument \Robject{n} for specifying the number of principal components to calculate.

%%%%%------- offline analysis --------%%%%%%

% <<fit-prcomp-offline>>=
% matter:::profile({
%   prcomp.out <- prcomp(data, n=2, retx=FALSE, center=FALSE, scale.=FALSE)
% })
%    # start (MB)   finish (MB) max used (MB) overhead (MB)    time (sec) 
%    #     84.200        84.600      1004.800       920.200       224.628 
%     # 84.200        84.600      1004.800       920.200       226.084 
%      # 84.200        84.600      1004.800       920.200       244.885 
%        # 84.200        84.600      1004.800       920.200       247.294 
%         # 84.400        84.800      1005.000       920.200       237.696 

% data.m <- data[]

% matter:::profile({
%   svd.out <- svd(data.m, nu=0, nv=2)
% })
%    # start (MB)   finish (MB) max used (MB) overhead (MB)    time (sec) 
%    #   1288.700      1288.700      3689.100      2400.400        68.158 
%     # 1288.700      1288.700      3689.100      2400.400        67.937 
%       # 1288.70       1288.70       3689.10       2400.40         67.53 
%      # 1288.700      1288.700      3689.100      2400.400        84.922 
%       # 1288.900      1288.900      3689.300      2400.400        73.907 

% matter:::profile({
%   irlba.out <- irlba(data.m, nu=0, nv=2, fastpath=FALSE)
% })
%    # start (MB)   finish (MB) max used (MB) overhead (MB)    time (sec) 
%    #   1300.900      1300.900      3546.400      2245.500        11.653
%     # 1288.700      1300.700      3540.100      2239.400        12.085  
%       # 1288.700      1300.700      3540.100      2239.400        11.855 
%      # 1288.700      1300.700      3540.100      2239.400        13.295 
%       # 1288.900      1300.900      3540.300      2239.400        12.294 

% require(irlba)
% require(bigalgebra)
% data.bm <- as.big.matrix(data.m)

% rm(data.m)
% gc()

% mult <- function(A, B, transpose=FALSE) {
%     if ( is.null(dim(B)) )
%       B <- as.matrix(B)
%     if ( is.null(dim(A)) )
%       A <- t(A)
%     if ( transpose ) {
%       cbind((t(B) %*% A)[])
%     } else {
%       cbind((A %*% B)[])
%     }
% }

% matter:::profile({
%   irlba.out <- irlba(data.bm, nu=0, nv=2, mult=mult)
% })
%    # start (MB)   finish (MB) max used (MB) overhead (MB)    time (sec) 
%    # 90.800       103.100      1926.100      1823.000        10.012 
%    # 90.600       103.100      1926.200      1823.100        10.183
%    # 90.600       103.100      1926.200      1823.100        10.737

% rm(data.bm)

% require(ffbase)
% require(bootSVD)   
% data.m <- data[]
% data.ff <- ff(vmode="double", dim=dim(data.m))
% data.ff[] <- data.m

% rm(data.m)
% gc()

% mult <- function(A, B, transpose=FALSE) {
%     if ( is.null(dim(B)) )
%       B <- as.matrix(B)
%     if ( is.null(dim(A)) )
%       A <- t(A)
%     if ( transpose ) {
%       cbind(ffmatrixmult(t(B), A)[])
%     } else {
%       cbind(ffmatrixmult(A, B)[])
%     }
% }

% matter:::profile({
%   irlba.out <- irlba(data.ff, nu=0, nv=2, mult=mult)
% })
%   # start (MB)   finish (MB) max used (MB) overhead (MB)    time (sec) 
%        # 89.200       102.200      1580.800      1478.600        83.821 
%        # 89.200       102.200      1021.000       918.800        80.844
%        # 89.600       102.300      1021.000       918.700        81.577
%        # 89.200       102.200      1021.000       918.800        91.519
%         # 89.200       102.200      1021.000       918.800        78.467
% @

%%%%%------- online analysis --------%%%%%%

<<fit-prcomp, eval=FALSE>>=
prcomp.out <- prcomp(data, n=2, center=FALSE, scale.=FALSE)
@

On a 2012 retina MacBook Pro with 2.6 GHz Intel CPU, 16 GB RAM, and 500 GB SSD, calculating the first two principal components takes roughly 2-4 minutes and uses an additional 337 MB of memory overhead. The max amount of memory used during the computation was only 484 MB, for the 1.2 GB dataset.

Now we plot the first two principal components.

<<plot-pc1>>=
plot(prcomp.out$rotation[,1], type='h', ylab="PC 1")
@

<<plot-pc2>>=
plot(prcomp.out$rotation[,2], type='h', ylab="PC 2")
@

\setkeys{Gin}{width=\textwidth}
\begin{figure}[h]
\centering
\begin{subfigure}{.2\textwidth}
  \centering
<<fig=TRUE, echo=FALSE, results=hide>>=
<<plot-var>>
@
\caption{\small Sample variance}
\label{fig:var}
\end{subfigure}
\begin{subfigure}{.2\textwidth}
  \centering
<<fig=TRUE, echo=FALSE, results=hide>>=
<<plot-pc1>>
@
\caption{\small PC1 loadings}
\label{fig:pc1}
\end{subfigure}
\begin{subfigure}{.2\textwidth}
  \centering
<<fig=TRUE, echo=FALSE, results=hide>>=
<<plot-pc2>>
@
\caption{\small PC2 loadings}
\label{fig:pc2}
\end{subfigure}
\caption{\small Principal components analysis for on-disk dataset.}
\end{figure}

As shown in the plots of the first and second principal components, the first PC shows that most of the variation in the data occurs in the first twenty variables, while the second PC distinguishes the first ten variables from the next ten variables.


%%%%%------- save data ----------%%%%%%

% <<save-data>>=
% save(bigglm.out, prcomp.out, file="~/Documents/Developer/Projects/matter/data/matter_ex.rda")
% @



\section{Design and implementation}

\Rpackage{matter} is designed with several goals in mind. Like the \textit{bigmemory} and \textit{ff} packages, it seeks to make statistical methods scalable to larger-than-memory datasets by utilizing data-on-disk. Unlike those packages, it seeks to make domain-specific file formats (such as Analyze 7.5 and imzML for MS imaging experiments) accessible from disk directly without additional file conversion. It seeks to have a minimal memory footprint, and require minimal developer effort to use, while maintaining computational efficiency wherever possible.

\subsection{S4 classes}

\Rpackage{matter} utilizes S4 classes to implement on-disk matrices in a way so that they can be seamlessly accessed as if they were ordinary \R{} matrices. These are the \Robject{atoms} class and the \Robject{matter} class. The \Robject{atoms} class is not exported to the user, who only interacts with the \Robject{matter} class and \Robject{matter} objects to create and manipulate on-disk matrices.

\subsubsection{\Robject{atoms}: contiguous sectors of data on disk}

By analogy to \R{}'s notion of ``atomic'' vectors, the \Robject{atoms} class uses the notion of contiguous ``atomic'' sectors of disk. Each ``atom'' in an \Robject{atoms} object gives the location of one block of contiguous data on disk, as denoted by a file path, an byte offset from the beginning of the file, a data type, and the number of data elements (i.e., the length) of the atom. An \Robject{atoms} object may consist of many atoms from multiple files and multiple locations on disk, while ultimately representing a single vector or row or column of a matrix.

Structure:
\begin{itemize}
\item \Robject{length}: the total cumulative length of all of the atoms
\item \Robject{source\_id}: the ID's of the file paths where each atom is located
\item \Robject{datamode}: the type of data (short, int, long, float, double) for each atom
\item \Robject{offset}: each atom's byte offset from the beginning of the file
\item \Robject{extent}: the length of each atom
\item \Robject{index\_offset}: the cumulative index of the first element of each atom
\item \Robject{index\_extent}: the cumulative one-past-the-end index of each atom
\end{itemize}

The \Robject{atoms} class has a C++ backend in the \Robject{Atoms} C++ class.

\subsubsection{\Robject{matter}: vectors and matrices stored on disk}

A \Robject{matter} object is made of one or more \Robject{atoms} objects, and represents a vector or matrix. It includes additional metadata such as dimensions and row names or column names.

Structure:
\begin{itemize}
\item \Robject{data}: one or more \Robject{atoms} objects
\item \Robject{datamode}: the type of data (integer, numeric) for the represented vector or matrix
\item \Robject{paths}: the paths to the files used by the \Robject{atoms} objects
\item \Robject{filemode}: should the files be open for read/write, or read-only?
\item \Robject{chunksize}: how large the chunk sizes should be for calculations that operate on chunks of the dataset
\item \Robject{length}: the total length of the dataset
\item \Robject{dim}: the extent of each dimension (for a matrix)
\item \Robject{names}: the names of the data elements (for a vector)
\item \Robject{dimnames}: the names of the dimensions (for a matrix)
\end{itemize}

A \Robject{matter\_vec} vector contains a single \Robject{atoms} object that represents all of the atoms of the vector. The \Robject{matter\_mat} matrix class has two subtypes for column-major (\Robject{matter\_matc}) and row-major (\Robject{matter\_matr}) matrices. A column-major \Robject{matter\_matc} matrix has one \Robject{atoms} object for each column, while a row-major \Robject{matter\_matr} matrix has one \Robject{atoms} object for each row.

The \Robject{matter} class has a C++ backend in the \Robject{Matter} C++ class.

\subsection{C++ classes}

\Rpackage{matter} utilizes a C++ backend to access the data on disk and transform it into the appropriate representation in R. Although these classes correspond to S4 classes in \R{}, and are responsible for most of the computational work, all of the required metadata is stored in the S4 classes in \R{}, and are simply read by the C++ classes. This means that \Rpackage{matter} never depends on external pointers, which makes it trivial to share \Robject{matter} vectors and \Robject{matter} matrices between \R{} sessions that have access to the same filesystem.

\subsubsection{\Robject{Atoms}: contiguous sectors of data on disk}

The \Robject{Atoms} C++ class is responsible for reading and writing the data on disk, based on its metadata. For computational efficiency, it tries to perform sequential reads over random read/writes whenever possible, while minimizing the total number of atomic read/writes on disk.

\subsubsection{\Robject{Matter}: vectors and matrices stored on disk}

The \Robject{Matter} C++ class is responsible for transforming the data read by the \Robject{Atoms} class into a format appropriate for R. This may include re-arranging contiguous data that has been read sequentially into a different order, either due to the inherent organization of the dataset, or as requested by the user in R.

\subsubsection{\Robject{MatterAccessor}: iterate over virtual disk objects}

The \Robject{MatterAccessor} C++ class acts similarly to an iterator, and allows buffered iteration over a \Robject{Matter} object. It can either iterate over the whole dataset (for both vectors and matrices), or over a single column for column-major matrices, or over a single row for row-major matrices.

A \Robject{MatterAccessor} object will load portions of the dataset (as many elements as the \Robject{chunksize} at once) into memory, and then free that portion of the data and load a new chunk, as necessary. This buffering is handled automatically by the class, and code can treat it as a regular iterator. This allows seamless and simple iteration over \Robject{Matter} objects while maintaining strict control over the memory footprint of the calculation.


\section{Session info}

<<results=tex, echo=FALSE>>=
toLatex(sessionInfo())
@

% \bibliographystyle{unsrt}
\bibliography{matter}

\end{document}
