# SCM Repository

[matrix] Diff of /pkg/R/bbCrosstab.R
 [matrix] / pkg / R / bbCrosstab.R

# Diff of /pkg/R/bbCrosstab.R

revision 347, Tue Nov 23 15:37:58 2004 UTC revision 348, Tue Nov 23 15:38:31 2004 UTC
# Line 1  Line 1
1  bbCrosstab <- function(flist) {  bCrosstab <- function(flist) {
2      ind <- function(i,j) ((i-1) * i)/2 + j # index in rowwise lower triangle      ind <- function(i,j) ((i-1) * i)/2 + j # index in rowwise lower triangle
3      flist <- lapply(as.list(flist), function (x) as(x, "factor")[drop = TRUE])      ## Coerce flist to a list of factors with no unused levels
4        flist <- lapply(as.list(flist), function (x) as.factor(x)[drop = TRUE])
5
6      nfac <- length(flist)      nfac <- length(flist)
7      if (nfac < 1) return(new("bbSparseSy", x = list(), uplo = 'L'))      if (nfac < 1) error("flist must be a non-empty list of factors")
8                                            # Check for consistent lengths
9        nobs <- length(flist[[1]])
10        if (any(lapply(flist, length) != nobs))
11            stop("all factors in flist must have the same length")
12        ones <- rep(1, nobs)
13        ## Order the factors by non-increasing length(levels(el))
14      nlev <- unlist(lapply(flist, function(x) length(levels(x))))      nlev <- unlist(lapply(flist, function(x) length(levels(x))))
15      ord <- rev(order(nlev))      ord <- rev(order(nlev))
16      if (any(ord != seq(along = ord))) {      if (any(ord != seq(along = ord))) {
17          nlev <- nlev[ord]          nlev <- nlev[ord]
18          flist <- flist[ord]          flist <- flist[ord]
19      }      }
20      nobs <- length(flist[[1]])                                          # establish output list and names
if (any(lapply(flist, length) != nobs))
stop("all factors in flist must have the same length")
ones <- rep(1, nobs)
21      lst <- vector("list", choose(nfac + 1, 2))      lst <- vector("list", choose(nfac + 1, 2))
22      if (is.null(nms <- names(flist))) nms <- paste("V", 1:nfac, sep = "")      if (is.null(nms <- names(flist))) nms <- paste("V", 1:nfac, sep = "")
23      nmat <- outer(nms, nms, FUN = "paste", sep = ":")      nmat <- outer(nms, nms, FUN = "paste", sep = ":")
# Line 20  Line 25
25      zb <- lapply(flist, function(x) as.integer(x) - 1:1) # zero-based indices      zb <- lapply(flist, function(x) as.integer(x) - 1:1) # zero-based indices
26      for (j in 1:nfac) {      for (j in 1:nfac) {
27          for (i in j:nfac) {          for (i in j:nfac) {
28              lst[[ind(i,j)]] <- as(as(new("tripletMatrix",              IND <- ind(i,j)
29                lst[[IND]] <- as(new("tripletMatrix",
30                                         i = zb[[i]], j = zb[[j]],                                         i = zb[[i]], j = zb[[j]],
31                                         x = ones, Dim = nlev[c(i,j)]),                                         x = ones, Dim = nlev[c(i,j)]),
32                                     "cscMatrix"),                               "cscMatrix")
"cscBlocked")
33              nms[ind(i,j)] <- nmat[i, j]              nms[ind(i,j)] <- nmat[i, j]
34          }          }
35      }      }
36      names(lst) <- nms      names(lst) <- nms
37      new("bbCrosstab", x = lst, uplo = "L")      list(flist = flist, ctab = lst)
38  }  }
39
40  setMethod("isNested", signature(object = "bbCrosstab"),  ctab2bblist <- function(ctab, nf, nc)
41            function(object, ...)  {
all(unlist(lapply(object@x, function(x) any(diff(x@p) > 1)))))

bbCrossNfac <- function(x) {            # determine the number of factors
x <- as(x, "bbCrosstab")
as.integer((-1 + sqrt(1 + 8*length(x@x)))/2)
}

bbCrossProject <- function(x) {         # Project the first column onto the others
42      ind <- function(i,j) ((i-1) * i)/2 + j # index in compressed lower triangle      ind <- function(i,j) ((i-1) * i)/2 + j # index in compressed lower triangle
43      x <- as(x, "bbCrosstab")      ans <- vector("list", length(ctab))
44      nf <- bbCrossNfac(x)      names(ans) <- names(ctab)
45      if (nf < 2) stop("Number of factors must be > 1 to project")      for (j in 1:nf) {
46      ans <- new("bbCrosstab", x = x@x[-(1:nf)], uplo = "L")          for (i in j:nf) {
47      for (i in 2:nf) {              IND <- ind(i, j)
48          ans@x[[ind(i-1,i-1)]] <-              ctij <- ctab[[IND]]
49              as(tcrossprod(as(x@x[[i]], "cscMatrix")), "cscBlocked")              ans[[IND]] <- new("cscBlocked", p = ctij@p, i = ctij@p,
50                                  x = array(0, dim = c(nc[i], nc[j], length(ctij@x))))
51            }
52      }      }
53      ans      ans
54  }  }
55
56  bbCrossPerms <- function(x) {  Lstruct <- function(bcr, nc = rep(1, nf)) {
57      ind <- function(i,j) ((i-1) * i)/2 + j      ind <- function(i,j) ((i-1) * i)/2 + j # index in compressed lower triangle
58      needsPerm <- function(x) any(diff(as(x, "cscBlocked")@p) != 1)      ctab <- bcr\$ctab
59      x <- as(x, "bbCrosstab")      fl <- bcr\$flist
60      nf <- bbCrossNfac(x)      nf <- length(fl)
61      perm <- vector("list", nf)      Linv <- vector("list", length(fl))
62  #    Lmat <- new("bbSparseTr", uplo = 'L', diag = 'U',      names(Linv) <- names(fl)
#                Linv = vector("list", nf), x = x@x)
Linv <- vector("list", nf)
for (j in 1:nf) {
db <- x@x[[ind(j,j)]]             # diagonal block
if (needsPerm(db)) {
res <- .Call("sscMatrix_ldl_symbolic",
as(as(db, "cscMatrix"), "sscMatrix"), TRUE)
} else {
ncol <- length(db@p) - 1
res <- list(rep(-(1:1), ncol), NULL, seq(ncol) - 1:1)
}
perm[[j]] <- res[[3]]
Linv[[j]] <- as(.Call("Parent_inverse", res[[1]], TRUE),
"cscBlocked")
}
#    Lmat\$Linv <- Linv
#    list(perm = perm, Lmat = Lmat)
list(perm = perm, Linv = Linv)
}
63
64        if (all(unlist(lapply(ctab, function(x) all(diff(x@p)== 1))))) { # nested
65            ZtZ <- ctab2bblist(ctab, nf, nc)
66            for (j in 1:nf) Linv[[j]] <- ZtZ[[ind(j, j)]]
67            return(list(flist = fl, ZtZ = ZtZ, Lmat = ZtZ, Linv = Linv))
68        }
69        ## non-nested case - here things get interesting
70        ct11 <- ctab[[1]]
71        Linv[[1]] <- new("cscBlocked", p = ct11@p, i = ct11@i,
72                         x = array(0, dim = c(nc[1], nc[1], length(ct11@x))))
73        for (i in 1:(nf - 1)) {
74            ip1 <- i + 1
75            res <- .Call("bCrosstab_project", ctab, i)
76            fac <- fl[[ip1]]
77            fl[[ip1]] <- factor(as.character(fac), levels = levels(fac)[1+res\$perm])
78            ctab <- .Call("bCrosstab_permute", res\$ctab, i, res\$perm)
79            rLi <- res\$Linv
80            Linv[[ip1]] <- new("cscBlocked", p = rLi@p, i = rLi@i,
81                               x = array(0, dim = c(nc[ip1], nc[ip1], length(rLi@x))))
82        }
83        list(flist = fl, ZtZ = ctab2bblist(bCrosstab(fl)\$ctab, nf, nc),
84             Lmat = ctab2bblist(ctab, nf, nc), Linv = Linv)
85        ## FIXME: You probably don't want to define Lmat from ctab because
86        ## it has off-diagonal elements generated from the inverses of the diagonals
87    }

Legend:
 Removed from v.347 changed lines Added in v.348