%\VignetteIndexEntry{gRbase-arrayops: Array operations in gRbase}
%\VignettePackage{gRbase}

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

\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}

\SweaveOpts{keep.source=T}

<<echo=FALSE,print=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)
##load_all("../gRbase")
options("width"=100)
options(useFancyQuotes="UTF-8")
##options(prompt=" ", "continue"="   ")
@ %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}


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

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.

Presence of named dimnames 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 


\subsection{Defining arrays}

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. 
Another way is to use \rr{newar()} from \grbase:
@ 
<<>>=
counts <- c(32, 53, 11, 50, 10, 25, 36, 66, 9, 34, 5, 29)
z3 <- newar( ~ Hair:Eye:Sex, levels = dn, value = counts) 
z4 <- newar(c("Hair", "Eye", "Sex"), levels=dn, values=counts)
@ %def 

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


Using \rr{newar()}, 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:
@ 
<<>>=
z6 <- newar(~Hair:Eye:Sex, levels=c(2, 3, 2), values = counts, normalize="first")
z6 %>% flat
@ %def 
The same can be achived with\rr{arnormalize()}
@ 
<<>>=
arnormalize(z5, "first") %>% flat
@ %def 


\section{Operations on single 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{Permuting an array}
\label{sec:permuting-an-array}

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

Alternative forms (that will also work for \rr{aperm()}):
@ 
<<results=hide>>=
arperm(hec, c("Eye", "Sex", "Hair")) %>% flat
arperm(hec, c(2,3,1)) %>% flat
@ %def 

Notice that abbreviation is allowed
@ 
<<results=hide>>=
arperm(hec, ~Ey:Se:Ha) %>% flat 
arperm(hec, c("Ey", "Se", "Ha")) %>% 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. Using standard \R\ code we can do:

@ 
<<>>=
hec[, 2:3, ]  %>% flat
is.array( hec[, 2:3, ] ) 
hec[1, , 1]
is.array( hec[1, , 1] )
@ %def 

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


\grbase\ provides a wrapper for this, for example:
@ 
<<results=hide>>=
arslice_prim(hec, slice=list(TRUE, 2:3, TRUE))  %>% flat
@ %def 


Using \rr{arslice()} from \grbase\ the following are equivalent
@ 
<<results=verbatim>>=
hec[, 2:3, ] %>% flat
arslice(hec, slice=list(c(2,3)), margin=2) %>% flat
@ %def 

Alternative forms:
@ 
<<results=hide>>=
arslice(hec, slice=list(Eye=2:3)) %>% flat
arslice(hec, slice=list(Eye=c("Blue","Hazel")))  %>% flat
@ %def 


The virtue of \rr{arslice} comes from the flexibility when
specifying the slice:\footnote{Currently names can not be abbreviated,
but that might be added later.}
@ 
<<>>=
arslice(hec, slice=list(Eye=c(2,3), Sex="Female")) %>% flat
arslice(hec, slice=list(Eye=c(2,3), Sex="Female"), drop=FALSE) %>% flat
@ %def 

If slicing leads to a one dimensional array, the output will by
default not be an array but a vector (without a dim attribute)
@ 
<<>>=
z <- arslice(hec, slice=list(Hair=1, Sex="Female")); z 
is.array( z )
@ %def 

The output can be forced to be an array with
@ 
<<>>=
z <- arslice(hec, slice=list(Hair=1, Sex="Female"), as.array=TRUE); z
is.array( z )
@ %def 




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

The $HE$--marginal table $T_{HE}$ is defined to be the table with
values
\begin{displaymath}
  T_{HE}(h,e) = \sum_s T_{HES}(h,e,s)
\end{displaymath}

With \grbase\ we can do\footnote{FIXME: Should allow for abbreviations in
  formula and character vector specifications.}:
@
<<>>=
he <- armarg(hec, ~Hair:Eye); he %>% flat       
@ %def

Alternative forms include
@ 
<<results=hide>>=
hs <- armarg(hec, c("Hair","Sex")); hs 
es <- armarg(hec, c(2,3)); es
@ %def 


\subsection{Inflating arrays}
\label{sec:inflating-arrays}

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_{HES}(h,e,s) = T_{HE}(h,e)
\end{displaymath}
so $\tilde T_{HES}(h,e,s)$ is constant as a function of $s$. 
With \grbase\ this is done with \rr{arexpand()}:
@ 
<<>>=
arexpand(he, list(Sex=c("Male", "Female"))) %>% flat
@ %def 
Notice that the added dimensions come ``at the end''. The following versions produce the same result:
@ 
<<results=hide>>=
arexpand(he, dimnames(hs)) %>% flat
arexpand(he, hs) %>% flat
@ %def 



\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}

With \grbase\ this is done with \rr{armult()}:
@ 
<<>>=
armult(he, hs) %>% flat
@ %def 


The quotient, sum and difference is defined similarly:
@ 
<<>>= 
ardiv(he, hs)  %>% flat 
aradd(he, hs)  %>% flat
arsubt(he, hs) %>% flat
@ %def 

Multiplication and addition of a list of arrays is accomplished with
\rr{arprod()} and \rr{arsum()}:
@ 
<<>>=
arprod( he, hs, es ) %>% flat
arsum( he, hs, es ) %>% flat
@ %def 


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

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 <- arperm(hec, 3:1)
arequal(hec, hec2)
@ %def 

A visual comparison of the entries of two arrays with the same
dimnames is much easier if the dimnames are in the same order.
For example, the following provides the fitted cell counts under a
specific log--linear model:
@ 
<<>>=
hec3 <- ardiv( armult( he, es ), armarg( hec, "Eye" ) )
hec3 %>% flat
@ %def 
Comparing these with the observed data is tricky because of the ordering:
@ 
<<>>=
hec %>% flat
@ %def 

The function \rr{aralign()} 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} 
@ 
<<>>=
aralign(hec3, hec)  %>% flat
@ %def 





\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 <- newar( ~rain, levels = lev, values = c(.2, .8) )
s_r <- newar( ~sprinkler:rain, levels = lev, values = c(.01,.99, .4, .6) )
w_sr <- newar( ~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{arprod()}:
ways:
@
<<>>=
joint <- arprod( 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)$ with \rr{ardist()}
@
<<>>=
rw <- armarg(joint, ~rain+wet)
ardist(rw, cond=~wet)
@ %def

@ 
<<eval=F>>=
## Alternative:
ardiv( rw, armarg(rw, ~wet))
## or
rw %a/% (rw %am% ~wet)
@ %def 


Alternative computation using 
@
<<>>=
x <- arslice_mult(rw, slice=list(wet="y")); x
p <- armarg(x, ~rain); p  ## Unnormalized
ardist( p )           ## Normalized
@ %def

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

We consider the $3$--way \code{lizard} data from \grbase:
@
<<>>=
data( lizard, package="gRbase" )
lizard
@ %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) armarg(indata, g))

    for (i in 1:4){
        for (j in seq_along(glist)){
            mf  <- armarg(fit, glist[[j]])
            ## adj <- ardiv( md[[ j ]], mf)
            ## fit <- armult( 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{nextCell()} }


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{nextCell()}.
@
<<>>=
nextCell(c(2,3,1), dim=c( 2, 3, 2 ))
@ %def




\subsection{\code{nextCellSlice()} 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{nextCellSlice()}.\footnote{FIXME: sliceset should be called margin.}
@ 
<<>>=
nextCellSlice(c(1,3,1), sliceset=2, dim=c( 2, 3, 2 ))
nextCellSlice(c(2,3,1), sliceset=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(slicecell=3, sliceset=2, dim=c( 2, 3, 2 ))
@ %def

To verify that we indeed get the right cells:

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


% 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



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

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

@ 
<<>>=
head( factGrid( 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 


% @
% <<>>=
% 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 





\end{document}




















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


% 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{nextCellSlice()}.
% @
% <<>>=
% nextCellSlice(c(2,1,1,2),  sliceset=2, dim2323)
% nextCellSlice(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{nextCell()}.

% @
% <<>>=
% nextCell(c(1,1,2,1), dim2222)
% nextCell(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{nextCellSlice()}.
% @
% <<>>=
% nextCellSlice(c(2,1,1,2),  sliceset=2, dim2323)
% nextCellSlice(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
