SCM

SCM Repository

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

Diff of /pkg/R/bbCrosstab.R

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge