%\VignetteEngine{knitr::knitr} 
%\VignetteIndexEntry{arrays: Array operations in gRbase}
%\VignettePackage{gRbase}

\documentclass[10pt]{article}
%\usepackage[T1]{fontenc}
%\usepackage[latin1]{inputenc}
%\usepackage{inputenx}
\usepackage{boxedminipage,color,a4wide,url}
\usepackage[T1]{fontenc}


\def\code#1{\texttt{#1}}
\def\R{\texttt{R}}
\def\pkg#1{\texttt{#1}}

\def\grain{\texttt{gRain}}
\def\grbase{\texttt{gRbase}}
\def\ptab{\code{ptab}}
\def\rr#1{\code{#1}[{\scriptsize gRbase}]}

\usepackage{fancyvrb}

\newlength{\fancyvrbtopsep}
\newlength{\fancyvrbpartopsep}
\makeatletter
\FV@AddToHook{\FV@ListParameterHook}{\topsep=\fancyvrbtopsep\partopsep=\fancyvrbpartopsep}
\makeatother

\setlength{\fancyvrbtopsep}{0pt}
\setlength{\fancyvrbpartopsep}{0pt}


\usepackage{etoolbox} 
\makeatletter 
\preto{\@verbatim}{\topsep=-10pt \partopsep=-10pt } 
\makeatother

<<include=FALSE>>=
library(knitr)
opts_chunk$set(
tidy=FALSE,
size="small"
)
@


<<echo=FALSE>>=
#require( gRbase )
prettyVersion <- packageDescription("gRbase")$Version
prettyDate <- format(Sys.Date())
@

\title{Array operations in the \grbase\ package}
\author{S{\o}ren H{\o}jsgaard}
\date{\pkg{gRbase} version \Sexpr{prettyVersion} as of \Sexpr{prettyDate}}


\begin{document}

%%\SweaveInput{Rmarkup.STY}

\definecolor{darkred}{rgb}{.7,0,0}
\definecolor{midnightblue}{rgb}{0.098,0.098,0.439}
%% 
%% \DefineVerbatimEnvironment{Sinput}{Verbatim}{
%%   fontfamily=tt,
%%   %%fontseries=b,
%%   %% xleftmargin=2em,
%%   formatcom={\color{midnightblue}}
%% }
%% \DefineVerbatimEnvironment{Soutput}{Verbatim}{
%%   fontfamily=tt,
%%   %%fontseries=b,
%%   %% xleftmargin=2em,
%%   formatcom={\color{darkred}}
%% }
%% \DefineVerbatimEnvironment{Scode}{Verbatim}{
%%   fontfamily=tt,
%%   %%fontseries=b,
%%   %% xleftmargin=2em,
%%   formatcom={\color{blue}}
%% }
%% 
%%\fvset{listparameters={\setlength{\topsep}{-2pt}}}
%%\renewenvironment{Schunk}{\linespread{.95}}{}

\maketitle

\parindent0pt\parskip5pt

\tableofcontents

<<echo=F>>=
library(gRbase)
options("width"=100, "digits"=4)
options(useFancyQuotes="UTF-8")
#chk = 'markup'
chk = 'hide'
@ %def


\section{Introduction}

This note describes some operations on arrays in \R. These operations
have been implemented to facilitate implementation of graphical models
and Bayesian networks in \R.

\section{Arrays/tables in \R}
\label{sec:arrays}

The documentation of \R\ states the following about arrays:

\begin{quote}
  \em An array in R can have one, two or more dimensions. It is simply
  a vector which is stored with additional attributes giving the
  dimensions (attribute "dim") and optionally names for those
  dimensions (attribute "dimnames").  A two-dimensional array is the
  same thing as a matrix.  One-dimensional arrays often look like
  vectors, but may be handled differently by some functions.
\end{quote}



%% <<>>=
%% apropos("^tab[A-Z]", ignore.case=FALSE)
%% @ %def

%% <<>>=
%% apropos("^table[A-Z]", ignore.case=FALSE)
%% @ %def

%% <<>>=
%% apropos("^ar_", ignore.case=FALSE)
%% @ %def


\subsection{Cross classified data - contingency tables}
\label{sec:new}

Arrays appear for example in connection with cross classified data. The array
\code{hec} below is an excerpt of the \code{HairEyeColor} array in \R:

<< >>= 
hec <- c(32, 53, 11, 50, 10, 25, 36, 66, 9, 34, 5, 29) 
dim(hec) <- c(2, 3, 2)
dimnames(hec) <- list(Hair = c("Black", "Brown"), 
                      Eye = c("Brown", "Blue", "Hazel"), 
                      Sex = c("Male", "Female"))
hec
@

Above, \code{hec} is an array because it has a \code{dim} attribute. Moreover,
\code{hec} also has a \code{dimnames} attribute naming the levels of each
dimension. Notice that each dimension is given a name.

Printing arrays can take up a lot of space.  Alternative views on an
array can be obtained with \code{ftable()} or by converting the array
to a dataframe with \code{as.data.frame.table()}. We shall do so in the following.

<<>>=
##flat <- function(x) {ftable(x, row.vars=1)}
flat <- function(x, n=4) {as.data.frame.table(x) %>% head(n)}
hec %>% flat
@ %def 

An array with named dimensions is in this package called a \emph{named array}.
The functionality described below relies heavily on arrays having named dimensions.
A check for an object being a named array is provided by
\rr{is.named.array()}
<<>>=
is.named.array(hec)
@ %def 


\subsection{Defining arrays}


Another way is to use \rr{tabNew()} from \grbase. This function is flexible wrt the input; for example:
<<>>=
dn <- list(Hair=c("Black", "Brown"), Eye=~Brown:Blue:Hazel, Sex=~Male:Female)
counts <- c(32, 53, 11, 50, 10, 25, 36, 66, 9, 34, 5, 29)
z3 <- tabNew(~Hair:Eye:Sex, levels=dn, value=counts) 
z4 <- tabNew(c("Hair", "Eye", "Sex"), levels=dn, values=counts)
@ %def 

Notice that the levels list (\code{dn} above) when used in \rr{tabNew()} 
is allowed to contain superfluous elements. Default
\code{dimnames} are generated with
<<>>=
z5 <- tabNew(~Hair:Eye:Sex, levels=c(2, 3, 2), values = counts)
dimnames(z5) %>% str
@ %def 


Using \rr{tabNew}, arrays can be normalized to sum to one in two ways:
1) Normalization can be over the first variable for \emph{each}
configuration of all other variables and 2) over all configurations. For
example:
<<>>=
z6 <- tabNew(~Hair:Eye:Sex, levels=c(2, 3, 2), values=counts, normalize="first")
z6 %>% flat
@ %def 


\section{Operations on arrays}
\label{sec:operations-arrays}

In the following we shall denote the dimnames 
(or variables) of the array \code{hec} by $H$, $E$ and $S$ and we let $(h,e,s)$
denote a configuration of these variables. The contingency table above
shall be denoted by $T_{HES}$ and we shall refer to the
$(h,e,s)$-entry of $T_{HES}$ as $T_{HES}(h,e,s)$. 

\subsection{Normalizing an array}
\label{sec:numarlizing-an-array}

Normalize an array with  \rr{tabNormalize()}
Entries of an array can be normalized to sum to one in two ways:
1) Normalization can be over the first variable for \emph{each}
configuration of all other variables and 2) over all configurations. For
example:
<<>>=
tabNormalize(z5, "first") %>% flat
@ %def 


\subsection{Subsetting an array -- slicing}
\label{sec:subsetting-an-array}

We can subset arrays (this will also be called ``slicing'') in
different ways. Notice that the result is not necessarily an
array. Slicing can be done using standard \R\ code or using \rr{tabSlice}.
The virtue of \rr{tabSlice} comes from the flexibility when
specifying the slice:


The following leads from the original $2\times 3 \times 2$
array to a $2 \times 2$
array by cutting away the \code{Sex=Male} and \code{Eye=Brown} slice of the array:
<<results=chk>>=
tabSlice(hec, slice=list(Eye=c("Blue", "Hazel"), Sex="Female"))
## Notice: levels can be written as numerics
## tabSlice(hec, slice=list(Eye=2:3, Sex="Female"))
@ 

We may also regard the result above as a $2 \times 2 \times 1$ array:
<<>>=
tabSlice(hec, slice=list(Eye=c("Blue", "Hazel"), Sex="Female"), drop=FALSE)
@ 

If slicing leads to a one dimensional array, the output will by
default not be an array but a vector (without a dim attribute). However, 
the result can be forced to be a 1-dimensional array:

<<results=chk>>=
## A vector:
t1 <- tabSlice(hec, slice=list(Hair=1, Sex="Female")); t1
## A 1-dimensional array:
t2 <- tabSlice(hec, slice=list(Hair=1, Sex="Female"), as.array=TRUE); t2 
## A higher dimensional array (in which some dimensions only have one level)
t3 <- tabSlice(hec, slice=list(Hair=1, Sex="Female"), drop=FALSE); t3
@ %def 

The difference between the last two forms can be clarified:
<<>>=
t2 %>% flat
t3 %>% flat
@ 



\subsection{Collapsing  and inflating arrays}
\label{sec:collapsing-arrays}

Collapsing: The $HE$--marginal array $T_{HE}$ of $T_{HES}$ is  the array with
values
\begin{displaymath}
  T_{HE}(h,e) = \sum_s T_{HES}(h,e,s)
\end{displaymath}
Inflating: The ``opposite'' operation is to extend an array. For example, we can
extend $T_{HE}$ to have a third dimension, e.g.\ \code{Sex}. That is
\begin{displaymath}
  \tilde T_{SHE}(s,h,e) = T_{HE}(h,e)
\end{displaymath}
so $\tilde T_{SHE}(s,h,e)$ is constant as a function of $s$. 

With \grbase\ we can collapse arrays with\footnote{FIXME: Should allow for abbreviations in
  formula and character vector specifications.}:
<<>>=
he <- tabMarg(hec, c("Hair", "Eye"))
he
@

<<results=chk>>=
## Alternatives
tabMarg(hec, ~Hair:Eye)
tabMarg(hec, c(1, 2))
hec %a_% ~Hair:Eye
@

Notice that collapsing is a projection in the sense that applying the
operation again does not change anything:
<<>>=
he1 <- tabMarg(hec, c("Hair", "Eye"))
he2 <- tabMarg(he1, c("Hair", "Eye"))
tabEqual(he1, he2)
@ 

Expand an array by adding additional dimensions with \rr{tabExpand()}:
<<>>=
extra.dim <- list(Sex=c("Male", "Female"))
tabExpand(he, extra.dim) 
@ %def 

<<results=chk>>=
## Alternatives
he %a^% extra.dim
@

Notice that expanding and collapsing brings us back to where we started:
<<>>=
(he %a^% extra.dim) %a_% c("Hair", "Eye")
@ 


\subsection{Permuting an array}
\label{sec:permuting-an-array}

A reorganization of the table can be made with \rr{tabPerm} (similar
to \code{aperm()}), but \rr{tabPerm} allows for a formula and for variable abbreviation:
<<>>=
tabPerm(hec, ~Eye:Sex:Hair) %>% flat 
@ %def 

Alternative forms (the first two also works for \code{aperm}):
<<results=chk>>=
tabPerm(hec, c("Eye", "Sex", "Hair"))
tabPerm(hec, c(2, 3, 1)) 
tabPerm(hec, ~Ey:Se:Ha) 
tabPerm(hec, c("Ey", "Se", "Ha"))
@ %def 




\subsection{Equality}
\label{sec:equality}

Two arrays are defined to be identical 1) if they have the same dimnames
and 2) if, possibly after a permutation, all values are identical (up to
a small numerical difference):

<<>>=
hec2 <- tabPerm(hec, 3:1)
tabEqual(hec, hec2)
@ %def 

<<results=chk>>= 
## Alternative
hec %a==% hec2
@

\subsection{Aligning}
\label{sec:aligning}

We can align one array according to the ordering of another:
<<results=chk>>= 
hec2 <- tabPerm(hec, 3:1)
tabAlign(hec2, hec)
@

<<result=chk>>=
## Alternative:
tabAlign(hec2, dimnames(hec))
@ 



%\section{Operations on two or more arrays}
%\label{sec:oper-two-arrays}

\subsection{Multiplication, addition etc:  $+$, $-$, $*$, $/$}
\label{sec:mult-addt-etc}

The product of two arrays $T_{HE}$ and $T_{HS}$ is defined to be the array
$\tilde T_{HES}$ with entries
\begin{displaymath}
  \tilde T_{HES}(h,e,s)= T_{HE}(h,e) + T_{HS}(h,s)
\end{displaymath}

The sum, difference and quotient is defined similarly: This is done
with \rr{tabProd()}, \rr{tabAdd()}, \rr{tabDiff()} and \rr{tabDiv()}:

<<>>=
hs <- tabMarg(hec, ~Hair:Eye)
tabMult(he, hs)
@ %def 

Available operations:
<<results=chk>>=
tabAdd(he, hs) 
tabSubt(he, hs)
tabMult(he, hs)
tabDiv(he, hs) 
tabDiv0(he, hs) ## Convention 0/0 = 0
@ %def 

Shortcuts:
<<results=chk>>= 
## Alternative
he %a+% hs
he %a-% hs
he %a*% hs
he %a/% hs
he %a/0% hs ## Convention 0/0 = 0
@

Multiplication and addition of (a list of) multiple arrays is
accomplished with \rr{tabProd()} and \rr{tabSum()} (much like
\rr{prod()} and \rr{sum()}):
<<>>=
es <- tabMarg(hec, ~Eye:Sex)
tabSum(he, hs, es)  
## tabSum(list(he, hs, es))
@ %def 


%% Lists of arrays are processed with
%% <<results=chk>>= 
%% tabListAdd(list(he, hs, es))
%% tabListMult(list(he, hs, es))
%% @

\subsection{An array as a probability density}
\label{sec:an-array-as}

If an array consists of non--negative numbers then it may be regarded as an
(unnormalized) discrete multivariate density. With this view, the following
examples should be self explanatory:
<< >>= 
tabDist(hec, marg=~Hair:Eye)
tabDist(hec, cond=~Sex) 
tabDist(hec, marg=~Hair, cond=~Sex) 
@

\subsection{Miscellaneous}
\label{sec:miscellaneous-1}

Multiply values in a slice by some number and all other values by
another number:
<< >>= 
tabSliceMult(es, list(Sex="Female"), val=10, comp=0)
@






%% 
%%For example, the following provides the fitted cell counts under a
%%specific log--linear model:
%%<<>>=
%%hec3 <- ar_div( ar_mult( he, es ), tabMarg( hec, "Eye" ) )
%%hec3 %>% flat
%%@ %def 
%%Comparing these with the observed data is tricky because of the ordering:
%%<<>>=
%%hec %>% flat
%%@ %def 
%%
%%The function \rr{ar\_align()} will align the first array to have the same
%%variable order as the second array which makes a visual comparison
%%easier:\footnote{FIXME: aralign should be modified so that the second
%%  argument can also be a list of dimnames} 
%%<<>>=
%%tabAlign(hec3, hec)  %>% flat
%%ar_perm(hec3, names(dimnames(hec)))  %>% flat
%%@ %def 
%%
%%<< >>= 
%%microbenchmark::microbenchmark(
%%                    tabAlign(hec3, hec),
%%                    ar_perm(hec3, names(dimnames(hec)))
%%                    )
%%@
%%


\section{Examples}
\label{sec:examples}


\subsection{A Bayesian network}
\label{sec:comp-with-arrays}

A classical example of a Bayesian network is the ``sprinkler
example'', see e.g.\
\url{http://en.wikipedia.org/wiki/Bayesian_network}:
\begin{quote}
  \em
  Suppose that there are two events which could cause grass to be wet:
  either the sprinkler is on or it is raining. Also, suppose that the
  rain has a direct effect on the use of the sprinkler (namely that
  when it rains, the sprinkler is usually not turned on). Then the
  situation can be modeled with a Bayesian network.
\end{quote}

We specify conditional probabilities $p(r)$, $p(s|r)$ and $p(w|s,r)$
as follows
(notice that the vertical conditioning bar ($|$) is replaced by the
horizontal underscore:

<<>>=
yn <- c("y","n")
lev <- list(rain=yn, sprinkler=yn, wet=yn)
r <- tabNew(~rain, levels=lev, values=c(.2, .8))
s_r <- tabNew(~sprinkler:rain, levels = lev, values = c(.01, .99, .4, .6))
w_sr <- tabNew( ~wet:sprinkler:rain, levels=lev, 
             values=c(.99, .01, .8, .2, .9, .1, 0, 1))
r 
s_r  %>% flat
w_sr %>% flat
@ %def 

The joint distribution $p(r,s,w)=p(r)p(s|r)p(w|s,r)$ can be obtained
with \rr{tabProd()}:
ways:
<<>>=
joint <- tabProd(r, s_r, w_sr); joint %>% flat
@ %def

What is the probability that it rains given that the grass is wet? We
find $p(r,w)=\sum_s p(r,s,w)$ and then $p(r|w)=p(r,w)/p(w)$. Can be done in various ways: with \rr{tabDist()}
<<>>=
tabDist(joint, marg=~rain, cond=~wet)
@ %def

<<results='hide'>>=
## Alternative:
rw <- tabMarg(joint, ~rain + wet)
tabDiv(rw, tabMarg(rw, ~wet))
## or
rw %a/% (rw %a_% ~wet)
@ %def 

<<>>=
## Alternative:
x <- tabSliceMult(rw, slice=list(wet="y")); x
tabDist(x, marg=~rain)
@ %def

\subsection{Iterative Proportional Scaling (IPS)}
\label{sec:ips}

We consider the $3$--way \code{lizard} data from \grbase:
<<>>=
data(lizard, package="gRbase")
lizard %>% flat
@ %def

Consider the two factor log--linear model for the \verb'lizard'
data. Under the model the expected counts have the form
\begin{displaymath}
  \log m(d,h,s)= a_1(d,h)+a_2(d,s)+a_3(h,s)
\end{displaymath}
If we let $n(d,h,s)$ denote the observed counts, the likelihood
equations are: Find $m(d,h,s)$ such that
\begin{displaymath}
  m(d,h)=n(d,h), \quad
  m(d,s)=n(d,s), \quad
  m(h,s)=n(h,s)
\end{displaymath}
where $m(d,h)=\sum_s m(d,h.s)$ etc. 
The updates are as follows: For the first term we have

\begin{displaymath}
  m(d,h,s) \leftarrow m(d,h,s) \frac{n(d,h)}{m(d,h)}
%  , \mbox{ where }
%  m(d,h) = \sum_s m(d,h,s)
\end{displaymath}
After iterating the updates will not change and we will have equality:
$  m(d,h,s) = m(d,h,s) \frac{n(d,h)}{m(d,h)}$ and summing over $s$
shows that the equation $m(d,h)=n(d,h)$ is satisfied. 

A rudimentary implementation of iterative proportional scaling for
log--linear models is straight forward:
<<>>=
myips <- function(indata, glist){
    fit   <- indata
    fit[] <-  1
    ## List of sufficient marginal tables
    md    <- lapply(glist, function(g) tabMarg(indata, g))

    for (i in 1:4){
        for (j in seq_along(glist)){
            mf  <- tabMarg(fit, glist[[j]])
            # adj <- tabDiv( md[[ j ]], mf)
            # fit <- tabMult( fit, adj )
            ## or
            adj <- md[[ j ]] %a/% mf
            fit <- fit %a*% adj
        }
    }
    pearson <- sum((fit - indata)^2 / fit)
    list(pearson=pearson, fit=fit)
}

glist <- list(c("species", "diam"),c("species", "height"),c("diam", "height"))

fm1 <- myips(lizard, glist)
fm1$pearson
fm1$fit %>% flat

fm2 <- loglin(lizard, glist, fit=T)
fm2$pearson
fm2$fit %>% flat
@ %def


\section{Some low level functions}
\label{sec:some-low-level}

For e.g.\ a $2\times 3 \times 2$ array, the entries are such that the first
variable varies fastest so the ordering of the cells are $(1,1,1)$,
$(2,1,1)$, $(1,2,1)$, $(2,2,1)$,$(1,3,1)$ and so on. To find the value
of such a cell, say,
$(j,k,l)$ in the array (which is really just a vector), the cell is
mapped into an entry of a vector. 

For example, cell $(2,3,1)$
(\verb|Hair=Brown|, \verb|Eye=Hazel|, \verb|Sex=Male|) must be mapped to
entry $4$ in
<<>>=
hec
c(hec)
@ %def 

For illustration we do:
<<>>=
cell2name <- function(cell, dimnames){
    unlist(lapply(1:length(cell), function(m) dimnames[[m]][cell[m]]))
}
cell2name(c(2,3,1), dimnames(hec))
@ %def 



\subsection{\code{cell2entry()}, \code{entry2cell()} and \code{next\_cell()} }


The map from a cell to the corresponding
entry is provided by \rr{cell2entry()}. The reverse operation, going
from an entry to a cell (which is much less needed) is provided by
\rr {entry2cell()}.

<<>>=
cell2entry(c(2,3,1), dim=c(2, 3, 2))
entry2cell(6, dim=c(2, 3, 2))
@ %def 


Given a cell, say $i=(2,3,1)$ in a $2\times 3\times 2$ array we often want to find the next cell in
the table following the convention that the first factor varies
fastest, that is $(1,1,2)$. This is provided by
\rr{next\_cell()}.
<<>>=
next_cell(c(2,3,1), dim=c(2, 3, 2))
@ %def


\subsection{\code{next\_cell\_slice()} and \code{slice2entry()}}
%\label{sec:x}

Given that we look at cells for which for which the index in dimension $2$ is at level $3$ (that is
\verb|Eye=Hazel|), i.e.\ cells of the form $(j,3,l)$. Given such a
cell, what is then the next cell that also satisfies this
constraint. This is provided by
\rr{next\_cell\_slice()}.\footnote{FIXME: sliceset should be called margin.}
<<>>=
next_cell_slice(c(1,3,1), slice_marg=2, dim=c( 2, 3, 2 ))
next_cell_slice(c(2,3,1), slice_marg=2, dim=c( 2, 3, 2 ))
@ %def 


Given that in dimension $2$ we look at level $3$. We want to find
entries for the cells of the form $(j,3,l)$.\footnote{FIXME:slicecell and
  sliceset should be renamed}
<<>>=
slice2entry(slice_cell=3, slice_marg=2, dim=c( 2, 3, 2 ))
@ %def

To verify that we indeed get the right cells:

<<>>=
r <- slice2entry(slice_cell=3, slice_marg=2, dim=c( 2, 3, 2 ))
lapply(lapply(r, entry2cell, c( 2, 3, 2 )),
       cell2name, dimnames(hec))
@ %def




\subsection{\code{fact\_grid()} -- Factorial grid}
\label{sec:factgrid}

Using the operations above we can obtain the combinations of the
factors as a matrix:

<<>>=
head( fact_grid( c(2, 3, 2) ), 6 )
@ %def 

A similar dataframe can also be obtained with the standard \R\
function \code{expand.grid} (but \code{factGrid} is faster)
<<>>=
head( expand.grid(list(1:2, 1:3, 1:2)), 6 )
@ %def 



\appendix

\section{More about slicing}
\label{sec:more-about-slicing}

Slicing using standard \R\ code can be done as follows:

<<>>=
hec[, 2:3, ]  %>% flat  ## A 2 x 2 x 2 array
hec[1, , 1]             ## A vector
hec[1, , 1, drop=FALSE] ## A 1 x 3 x 1 array
@ %def 

Programmatically we can do the above as
<<results=chk>>=
do.call("[", c(list(hec), list(TRUE, 2:3, TRUE)))  %>% flat
do.call("[", c(list(hec), list(1, TRUE, 1))) 
do.call("[", c(list(hec), list(1, TRUE, 1), drop=FALSE)) 
@ %def 

\grbase\ provides two alterntives for each of these three cases above:
<<results=chk>>=
tabSlicePrim(hec, slice=list(TRUE, 2:3, TRUE))  %>% flat
tabSlice(hec, slice=list(c(2, 3)), margin=2) %>% flat

tabSlicePrim(hec, slice=list(1, TRUE, 1))  
tabSlice(hec, slice=list(1, 1), margin=c(1, 3)) 

tabSlicePrim(hec, slice=list(1, TRUE, 1), drop=FALSE)  
tabSlice(hec, slice=list(1, 1), margin=c(1, 3), drop=FALSE) 
@ %def 

\end{document}




















%%%
%%%
%%% CUTOUT BELOW HERE
%%%
%%%


% Given $A\subset \Delta$ and a cell $i_A \in I_A$. This cell defines a
% slice of the original array, namely the cells
% $I(i_A)=\{j\in I|j_A = i_A\}$.
% We often want to find the entries in
% $x$ for the cells $I(i_A)$. This is provided by
% \code{slice2entry()}. For example, we may want the entries for
% the cells $(*,1,2,*)$ or $(2,2,*,*)$:

% @
% <<>>=
% r1<-slice2entry(slicecell=c(1,2), sliceset=c(2,3), dim2222); r1
% @ %def

% To verify that we indeed get the right cells:

% @
% <<>>=
% do.call(rbind, lapply(r1, entry2cell, dim2222))
% @ %def



% \subsection{\code{permuteCellEntries()}}

% %% SOMETHING WRONG WITH THE CODE; BUT WE DONT REALLY NEED IT

% For e.g.\ a $2\times 3 \times 2$ array, the entries $1,2,3, \dots, 12$
% correspond to the cells
% $(1,1,1)$,
% $(2,1,1)$, $(1,2,1)$, $(2,2,1)$,$(1,3,1)$, $\dots, (2,3,2)$. 

% Suppose dimensions are permuted as $(1,2,3) \rightarrow (2,3,1)$. Then
% the entries become:
% @ 
% <<>>=
% p <- permuteCellEntries(perm=c(2,1,3), dim=c( 2, 3, 2 ))
% p
% @ %def 

% @ 
% <<>>=
% hec[ p ]
% @ %def 




% In a $2\times 3$ table, entries $1,\dots,6$ correspond to combinations
% $(1,1),(2,1),(1,2),(2,2),(1,3),(2,3)$. If we permute the table to a $3
% \times 2$ table the entries become as follows:

% @
% <<>>=
% p<-permuteCellEntries(perm=c(2,1), dim=c(2,3)); p
% @ %def

% So for example,

% @
% <<>>=
% (A <- array(11:16, dim=c(2,3)))
% Ap <- A[p]
% dim(Ap) <- c(3,2)
% Ap
% @ %def

% This corresponds to

% @
% <<>>=
% aperm(A, c(2,1))
% @ %def

% @
% <<>>=
% ff <- factGrid(dim2222)
% head(ff, 4)
% tail(ff, 4)
% @ %def

% This is the same as (but faster)

% @
% <<>>=
% aa <- expand.grid(list(1:2,1:2,1:2,1:2))
% head(aa, 4)
% @ %def

% There is a slice version as well:
% @
% <<>>=
% factGrid(dim2222, slicecell=c(1,2), sliceset=c(2,3))
% @ %def

% @ 
% <<>>=
% if (require(microbenchmark)){
%     microbenchmark(as.data.frame(factGrid(c( 2, 3, 2 ))),
%                    expand.grid(list(1:2, 1:3, 1:2)))
% }
% @ %def 






% As an example we take the following:

% @
% <<>>=
% dim2222 <- c(2,2,2,2)
% dim2323 <- c(2,3,2,3)
% @ %def

% Given $A\subset \Delta$ and a cell $i_A \in I_A$ consider the cells
% $I(i_A)=\{j\in I|j_A = i_A\}$. For example, the cells satisfying that
% dimension $2$ is at level $1$. Given such a cell, say $(2,1,1,2)$ we
% often want to find the next cell also satisfying this constraint
% (again following the convention that the first factor varies fastest),
% that is $(1,1,2,2)$. This is provided by \code{next_cell_slice()}.
% @
% <<>>=
% next_cell_slice(c(2,1,1,2),  sliceset=2, dim2323)
% next_cell_slice(c(1,3,2,1),  sliceset=c(2,3), dim2323)
% @ %def



% Given a cell, say $i=(1,1,2,1)$ we often want to find the next cell in
% the table following the convention that the first factor varies
% fastest, that is $(2,1,2,1)$. This is provided by
% \code{next_cell()}.

% @
% <<>>=
% next_cell(c(1,1,2,1), dim2222)
% next_cell(c(2,2,2,1), dim2222)
% @ %def


% Given $A\subset \Delta$ and a cell $i_A \in I_A$ consider the cells
% $I(i_A)=\{j\in I|j_A = i_A\}$. For example, the cells satisfying that
% dimension $2$ is at level $1$. Given such a cell, say $(2,1,1,2)$ we
% often want to find the next cell also satisfying this constraint
% (again following the convention that the first factor varies fastest),
% that is $(1,1,2,2)$. This is provided by \code{next_cell_slice()}.
% @
% <<>>=
% next_cell_slice(c(2,1,1,2),  sliceset=2, dim2323)
% next_cell_slice(c(1,3,2,1),  sliceset=c(2,3), dim2323)
% @ %def




% @
% <<>>=
% cell2entry(c(1,1,1,1), dim2222)
% entry2cell(1, dim2222)
% cell2entry(c(2,1,2,1), dim2222)
% entry2cell(6, dim2222)
% @ %def





% @
% <<eval=F,echo=F>>=
% cell2entry2(c(2,1,2,1), dim)
% cell2entryR(c(2,1,2,1), dim)
% @ %def


%% @
%% <<>>=
%% ## because lizard is a vector we can do
%% lizard[1:2]
%% is.array( lizard[1:2] )
%% ## because lizard is an array we can do
%% lizard[,1,1]
%% is.array( lizard[,1,] )
%% @ %def


% @
% <<>>=
% data(lizard, package="gRbase")
% dimnames(lizard)
% @ %def

% @ 
% <<>>=
% lizard[, 1, ]
% is.array( lizard[, 1, ] )
% do.call("[", c(list(lizard), list(TRUE, 1, TRUE))) # programmatically
% tabSlice(lizard, slice=c(1), margin=c(2))          # using gRbase
% @ %def 

% If the result is a $1$--dimensional array, then dimensions are dropped
% by default unless one puts \code{drop=FALSE} in the appropriate places:
% @ 
% <<results=hide>>=
% lizard[, 1, 2]
% is.array( lizard[, 1, 2] )
% lizard[, 1, 2, drop=FALSE]
% is.array( lizard[, 1, 2, drop=FALSE] )
% do.call("[", c(list(lizard), list(TRUE, 1, 2), drop=FALSE))
% tabSlice(lizard, slice=c(1, 2), margin=c(2, 3), drop=FALSE)
% @ %def 



% \subsection{Slice of a table}
% \label{sec:slice-table}

% A slice of a table is obtained with \code{tabSlice}:

% @
% <<>>=
% tabSlice(lizard, slice=list(species="anoli"))
% tabSlice(lizard, slice=list(species="anoli", diam="<=4")) ## Dims are preserved...
% @ %def

% @ 
% <<>>=
% tabSlice(lizard, slice=list(species=c("anoli","dist"))) # FIXME: This fails.
% lizard[,,1:2]
% @ %def 


% \section{Operations on tables}
% \label{sec:operations-tables}


% Consider again the \verb'lizard' data In the following we shall denote
% the dimnames (or variables) by $D$, $H$ and $S$ and we let $(d,h,s)$
% denote a configuration of these variables. The contingency table above
% shall be denoted by $T_{DHS}$ and we shall refer to the
% $(d,h,s)$-entry as $T_{DHS}(d,h,s)$. 

% @
% <<>>=
% T.DHS <- lizard
% @ %def

% \subsection{Collapsing arrays onto marginals}
% \label{sec:marginal-tables}

% The $DS$--marginal table $T_{DS}$ is defined to be the table with
% values
% \begin{displaymath}
%   T_{DS}(d,s) = \sum_h T_{DHS}(d,h,s)
% \end{displaymath}

% @
% <<>>=
% T.DS <- tabMarg(lizard, ~diam+species); T.DS
% ## Alternative forms
% T.DS <- tabMarg(lizard, c("diam","species"))
% T.DS <- tabMarg(lizard, c(1,3))
% @ %def

% @
% <<echo=F,results=hide>>=
% T.DS <- tableMargin(lizard, ~diam+species); T.DS
% T.HS <- tableMargin(lizard, ~height+species); T.HS
% @ %def


% @ %def 
% @ 
% <<>>=
% a1 <- 1:12
% dim(a1) <- c( 2, 3, 2 )
% dimnames(a1) <- list(a=c("a1","a2"), b=c("b1","b2","b3"), c=c("c1","c2"))
% a1 %>% flat
% a2 <- array(1:12, dim=c( 2, 3, 2 ), 
%             dimnames=list(a=c("a1","a2"), b=c("b1","b2","b3"), c=c("c1","c2")))
% a2 %>% flat

% dim.names <- list(a=c("a1","a2"), b=c("b1","b2","b3"), c=c("c1","c2"), d=c("d1","d2"))
% tab(~a:b:c, dim.names, data=1:12) %>% flat
% ## tab(c("a","b","c"), dim.names, data=1:12) ## same thing
% tab(~a:b:c, c( 2, 3, 2 ), data=1:12) %>% flatle(row.vars=1)
% @ %def 


%% An alternative is to use 
%% \code{parray()} from \grbase. For example
%% @
%% <<>>=
%% yn <- c("y","n")
%% T.AB <- array(c(5,95,1,99), dim=c(2,2), dimnames=list("A"=yn, "B"=yn))
%% T.AB <- parray(c("A","B"), levels=list(yn, yn), values=c(5,95,1,99))
%% @ %def

%% Using \code{parray()}, arrays can be normalized in two ways:
%% Normalization can be over the first variable for \emph{each}
%% configuration of all other variables or over all configurations. For
%% example:

%% @
%% <<print=T>>=
%% T.AB <- parray(c("A","B"), levels=list(yn, yn), values=c(5,95,1,99),
%%                normalize="first")
%% T.AB <- parray(c("A","B"), levels=list(yn, yn), values=c(5,95,1,99),
%%                normalize="all")
%% @ %def

% Using \code{tab()}, arrays can be normalized in two ways:
% Normalization can be over the first variable for \emph{each}
% configuration of all other variables or over all configurations. For
% example:

% @ 
% <<>>=
% yn <- c("y","n")
% lev <- list("a"=yn, "b"=yn)
% tab(~a:b, levels=lev, data=c(5,95,1,99), normalize="first")
% tab(~a:b, levels=lev, data=c(5,95,1,99), normalize="all")
% @ %def 





% \subsection{Operations on two tables: $+$, $-$, $*$, $/$}
% \label{sec:oper-two-tabl}

% The product of two tables, e.g. $T_{DS}$ and $T_{HS}$ is defined to be
% the table $\tilde T_{DHS}$ with entries
% \begin{displaymath}
%   \tilde T_{DHS}(d,h,s)= T_{DS}(d,s) T_{HS}(h,s)
% \end{displaymath}

% In \R:
% @
% <<>>=
% t.ds <- tabMarg(lizard, ~diam:species)
% t.hs <- tabMarg(lizard, ~height:species)
% tabMult(t.ds, t.hs) %>% flat
% @ %def

% %% T.HS <- tabMarg(lizard, ~height+species)
% %% T.DHS.mult = tabMult( T.DS, T.HS )
% %% T.DHS.mult %>% flat


% %% @
% %% <<echo=F,results=hide>>=
% %% microbenchmark::microbenchmark(
% %%     tabMult__( T.DS, T.HS ), arrayOp( T.DS, T.HS ) )
% %% @ %def


% The quotient, sum and difference is defined similarly:
% %% @
% %% <<>>=
% %% T.DHS.div  = tabDiv( T.DS, T.HS )
% %% T.DHS.add  = tabAdd( T.DS, T.HS )
% %% T.DHS.subt = tabSubt( T.DS, T.HS )
% %% @ %def



% @ 
% <<>>=
% tabDiv(t.ds, t.hs) %>% flat
% tabAdd(t.ds, t.hs) %>% flat
% tabSubt(t.ds, t.hs) %>% flat
% @ %def 


%% @
%% <<>>=
%% flatle( tablePerm(T.DHS, c("species","height","diam")) )
%% @ %def


% \subsection{Expanding arrays to new dimensions}
% \label{sec:expanding-arrays-new}

% The ``opposite'' operation is to extend an array. For example, we can
% extend $T_{DS}$ to have a third dimension, e.g.\ \code{height}. That is
% \begin{displaymath}
%   T_{DHS}(d,h,s) = T_{DS}(d,s)
% \end{displaymath}
% so $T_{DHS}(d,h,s)$ is constant as a function of $h$. 

% @ 
% <<>>=
% tabExtend <- arrayExtendDomain  ## FIXME
% tabExtend(t.ds, list(height=c(">4.75", "<=4.75"))) %>% flat
% @ %def 



% @ 
% <<>>=
% tabExpand(t.ds, t.hs) %>% flat
% @ %def 




% Consider this way of ``blowing up'' an array with extra dimensions.
% @
% <<>>=
% T.HSD2 <- tabExpand(T.DS, T.HS); T.HSD2
% names(dimnames(T.HSD2))
% @ %def

% Here \verb'T.HSD2' is a $3$--way table with the same variable names as
% the union of the variable names in \verb'T.DS' and \verb'T.HS'. Those
% variables in those variables in \verb'T.HS' vary fastest. Lastly, if
% we regards \verb'T.HSD2' as a function of $(h,s,d)$ we see that
% \verb'T.HSD2' is constant as a function of $s$.


% \subsection{Cross classified data - contingency tables}

% Arrays appear for example in connection with cross classified data:
% Consider the \code{lizard} data in \grbase:
% @
% <<>>=
% data( lizard, package="gRbase" )
% lizard
% @ %def

% Data is a contingency table; a cross classified table of counts. In
% \R\ lingo, data is a \code{table} object, but it is also an array because
% it has a \code{dim} attribute:
% @ 
% <<>>=
% class( lizard )
% is.array( lizard )
% dim( lizard )
% @ %def 

% Alternative views of data are:
% @ 
% <<>>=
% flatle(lizard, row.vars=1)
% flat <- function(x) flatle(x, row.vars=1)
% lizard %>% flat # because gRbase imports the pipe %>% from magrittr
% @ %def 


%% @
%% <<>>=
%% x0 <- x1 <- x2 <- x3 <- 1:8
%% c(is.atomic(x0), is.vector(x0), is.matrix(x0), is.array(x0))

%% ## 1-dimensional array
%% dim(x1) <- 8
%% x1
%% c(is.atomic(x1), is.vector(x1), is.matrix(x1), is.array(x1))

%% ## 2-dimensional array (matrix)
%% dim(x2) <- c(2,4)
%% x2
%% c(is.atomic(x2), is.vector(x2), is.matrix(x2), is.array(x2))

%% ## 3-dimensional array
%% dim(x3) <- c(2,2,2) 
%% x3
%% c(is.atomic(x0), is.vector(x0), is.matrix(x0), is.array(x0))
%% @ %def

%% Notice that arrays do not need a \code{dimnames} attribute. However,
%% for some of the operations described in the following, \code{dimnames}
%% are essential. One may set dimnames with
%% @ 
%% <<>>=
%% dimnames(x3) <- list(a=c("a1","a2"), b=c("y","n"), c=c(1,2))
%% @ %def 



%% Hence the defining characterstic of an array is that it is a vector
%% with a dim attribute. For example
%% @
%% <<>>=
%% ## 1-dimensional array
%% x1 <- 1:8
%% dim(x1) <- 8
%% x1
%% c(is.array(x1), is.matrix(x1))

%% ## 2-dimensional array (matrix)
%% x2 <- 1:8
%% dim(x2) <- c(2,4)
%% x2
%% c(is.array(x2), is.matrix(x2))

%% ## 3-dimensional array
%% x3 <- 1:8
%% dim(x3) <- c(2,2,2)
%% x3
%% c(is.array(x3), is.matrix(x3))
%% @ %def



%%
%%Arrays appear for example in connection with cross classified
%%data. For example
%%<<>>=
%%HairEyeColor
%%@ %def 
%%
%%Data is a contingency table; a cross classified table of counts. In
%%\R\ lingo, data is a \code{table} object, but it is also an array because
%%it has a \code{dim} attribute:
%%<<>>=
%%class( HairEyeColor )
%%is.array( HairEyeColor )
%%dim( HairEyeColor )
%%@ %def 
%%
%%The array also has a
%%\code{dimnames} attribute and the list of \verb'dimnames' has names:
%%<<>>=
%%dimnames( HairEyeColor )
%%@ %def 
%%
%%
%%Notice from the output above that the first variable (\code{Hair})
%%varies fastest. The \verb'dimnames' attributes are important for many of the functions from
%%\grbase\ described in the following sections.
%%
%%\subsection{Named arrays}
%%\label{sec:named-arrays}
%%
%%An array with named dimnames is in this package called a \emph{named array}; this can be checked with
%%\rr{is.named.array()}
%%<<>>=
%%is.named.array( HairEyeColor )
%%@ %def 
%%
%%
%%
%%To limit output we shall only consider two hair colours and three eye colours.
%%<<>>=
%%hec <- do.call("[", c(list(HairEyeColor), list(1:2, 1:3, TRUE), drop=FALSE))
%%hec <- HairEyeColor[1:2, 1:3, ]
%%hec
%%@ %def 
%%
%%A more compact view of data can be achieved with \code{ftable()}.
%%Since \grbase\ imports the pipe operator \verb'%>%' from
%%the \pkg{magrittr} package we will in this note do:
%%<<>>=
%%flat <- function(x) {ftable(x, row.vars=1)}
%%hec %>% flat
%%@ %def 
%%

%% Arrays can be defined in different ways using standard \R\ code:
%% <<>>=
%% z1 <- c(32, 53, 11, 50, 10, 25, 36, 66, 9, 34, 5, 29)
%% di <- c(2, 3, 2)
%% dn <- list(Hair = c("Black", "Brown"), 
%%            Eye = c("Brown", "Blue", "Hazel"), 
%%            Sex = c("Male", "Female"))
%% dim(z1) <- di
%% dimnames(z1) <- dn
%% z2 <- array(c(32, 53, 11, 50, 10, 25, 36, 66, 9, 34, 5, 29),
%%             dim=di, dimnames=dn)
%% @ %def 
%% where the \code{dimnames} part in both cases is optional. 
