################################################################################
# 
#  mcatall:  Computation of Multiple CA (indicator matrix)
#
#  Input: D - Original data matrix of the form "objects x variables"
#
#  A modified MCA algorithm is applied to a Burt matrix
#  of the form "variable categories x variable categories". 
#  The algorithm results in the standard and principal coordinates of the row and column profiles 
#  of the corresponding Indicator Matrix.
#
#  The results are the same with mcatrad.m
#  
#
#  Output:     
#       Sr - Row standard coordinates
#	Sc - Column standard coordinates
#	Pr - Row principal coordinates 
#	Pc - Column principal coordinates
#
# Example:
#
# load talldata.txt, a random 267,264 x 6 dataset - place talldata.txt file under C:/ folder
# D <- read.table("C:/talldata.txt", header = TRUE, colClasses = "factor")
#
################################################################################

mcatall <- function(D)
{

# size of the original data matrix
I       <- dim(D)[1]
Q       <- dim(D)[2]

# indicator matrix

lev.n  <- unlist(lapply(D, nlevels))
n      <- cumsum(lev.n)
J.t    <- sum(lev.n)
Q.t    <- dim(D)[2]
Z      <- matrix(0, nrow = I, ncol = J.t)
newD <- lapply(D, as.numeric)
offset <- (c(0, n[-length(n)]))
for (i in 1:Q.t)
   Z[1:I + (I * (offset[i] + newD[[i]] - 1))] <- 1
fn     <- rep(names(D), unlist(lapply(D, nlevels)))
ln     <- unlist(lapply(D, levels))

dimnames(Z)[[2]] <- paste(fn, ln, sep = "")
dimnames(Z)[[1]] <- as.character(1:I)

B <- t(Z) %*% Z
# number of categories
J <- dim(B)[1]
# retain all non trivial dimensions
nd.max <- min(J-Q, I-1)

Fb     <- B / sum(B)
cmb    <- apply(Fb, 2, sum)
ePb    <- cmb %*% t(cmb)
Rb     <- (Fb - ePb) / sqrt(ePb)
decb  <- svd(Rb)
lamb   <- decb$d[1:nd.max]

# column standard / principal coordinates
Scz <- decb$u[,1:nd.max] / sqrt(cmb)
Pcz <- Scz[,1:nd.max] %*% diag(sqrt(lamb))

# row standard / principal coordinates
Prz <- Z%*%Scz*Q^(-1)
Srz <- Prz %*% diag(sqrt(lamb)^(-1))

# In case you want to fix signs
#Sr.sign <- sign(Srz[1,])
#Srz.sign     <- sign(Sr[1,])
#s.tr     <- rep(1, nd.max)
#s.tr[Sr.sign[1:nd.max] != Srz.sign] <- -1

#Sc.sign <- sign(Scz[1,])
#Scz.sign     <- sign(Sc[1,])
#sc.tr     <- rep(1, nd.max)
#sc.tr[Sc.sign[1:nd.max] != Scz.sign] <- -1

#Srz <- Srz %*%diag(s.tr)
#Prz <- Prz %*%diag(s.tr)

#Scz <- Scz %*%diag(sc.tr)
#Pcz <- Pcz %*%diag(sc.tr)

}