# SCM Repository

[matrix] View of /pkg/R/bCrosstab.R
 [matrix] / pkg / R / bCrosstab.R

# View of /pkg/R/bCrosstab.R

Mon Feb 7 09:25:50 2005 UTC (14 years, 5 months ago) by maechler
File size: 3518 byte(s)
`stop() for error()`
```bCrosstab <- function(flist, reorder = TRUE) {
ind <- function(i,j) ((i-1) * i)/2 + j # index in rowwise lower triangle
## Coerce flist to a list of factors with no unused levels
flist <- lapply(as.list(flist), function (x) as.factor(x)[drop = TRUE])

nfac <- length(flist)
if (nfac < 1) stop("flist must be a non-empty list of factors")
# Check for consistent lengths
nobs <- length(flist[[1]])
if (any(lapply(flist, length) != nobs))
stop("all factors in flist must have the same length")
ones <- rep(1, nobs)
nlev <- unlist(lapply(flist, function(x) length(levels(x))))
ord <- rev(order(nlev))
if (reorder && any(ord != seq(along = ord))) {
nlev <- nlev[ord]
flist <- flist[ord]
}
# establish output list and names
lst <- vector("list", choose(nfac + 1, 2))
if (is.null(nms <- names(flist))) nms <- paste("V", 1:nfac, sep = "")
nmat <- outer(nms, nms, FUN = "paste", sep = ":")
nms <- vector("character", length(lst))
zb <- lapply(flist, function(x) as.integer(x) - 1:1) # zero-based indices
for (j in 1:nfac) {
for (i in j:nfac) {
IND <- ind(i,j)
lst[[IND]] <- as(new("dgTMatrix",
i = zb[[i]], j = zb[[j]],
x = ones, Dim = nlev[c(i,j)]),
ifelse(j != i,"dgCMatrix","dsCMatrix"))
nms[ind(i,j)] <- nmat[i, j]
}
}
names(lst) <- nms
list(flist = flist, ctab = lst)
}

ctab2bblist <- function(ctab, nf, nc)
{
ind <- function(i,j) ((i-1) * i)/2 + j # index in compressed lower triangle
ans <- vector("list", length(ctab))
names(ans) <- names(ctab)
for (j in 1:nf) {
for (i in j:nf) {
IND <- ind(i, j)
ctij <- ctab[[IND]]
ans[[IND]] <- new("dgBCMatrix", p = ctij@p, i = ctij@p,
x = array(0, dim = c(nc[i], nc[j], length(ctij@x))))
}
}
ans
}

Lstruct <- function(bcr, nc = rep(1, nf)) {
ind <- function(i,j) ((i-1) * i)/2 + j # index in compressed lower triangle
ctab <- bcr\$ctab
fl <- bcr\$flist
nf <- length(fl)
Linv <- vector("list", length(fl))
names(Linv) <- names(fl)

if (all(unlist(lapply(ctab, function(x) all(diff(x@p)== 1))))) { # nested
ZtZ <- ctab2bblist(ctab, nf, nc)
for (j in 1:nf) Linv[[j]] <- ZtZ[[ind(j, j)]]
return(list(flist = fl, ZtZ = ZtZ, Lmat = ZtZ, Linv = Linv))
}
## non-nested case - here things get interesting
ct11 <- ctab[[1]]
Linv[[1]] <- new("dgBCMatrix", p = ct11@p, i = ct11@i,
x = array(0, dim = c(nc[1], nc[1], length(ct11@x))))
for (i in 1:(nf - 1)) {
ip1 <- i + 1
res <- .Call("bCrosstab_project", ctab, i)
fac <- fl[[ip1]]
fl[[ip1]] <- factor(as.character(fac), levels = levels(fac)[1+res\$perm])
ctab <- .Call("bCrosstab_permute", res\$ctab, i, res\$perm)
rLi <- res\$Linv
Linv[[ip1]] <- new("dgBCMatrix", p = rLi@p, i = rLi@i,
x = array(0, dim = c(nc[ip1], nc[ip1], length(rLi@x))))
}
list(flist = fl, ZtZ = ctab2bblist(bCrosstab(fl)\$ctab, nf, nc),
Lmat = ctab2bblist(ctab, nf, nc), Linv = Linv)
## FIXME: You probably don't want to define Lmat from ctab because
## it has off-diagonal elements generated from the inverses of the diagonals
}
```