SCM

SCM Repository

[matrix] Diff of /pkg/Matrix/R/diagMatrix.R
ViewVC logotype

Diff of /pkg/Matrix/R/diagMatrix.R

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

revision 2197, Mon Jun 2 14:34:42 2008 UTC revision 2226, Mon Jul 21 17:15:17 2008 UTC
# Line 32  Line 32 
32      }      }
33  }  }
34    
35  ## Pkg 'spdep' had (relatively slow) versions of this as_dsCMatrix_I()  .sparseDiagonal <- function(n, x = rep.int(1,n), uplo = "U", shape = "t") {
 .symDiagonal <- function(n, x = rep.int(1,n), uplo = "U") {  
36      stopifnot(n == (n. <- as.integer(n)), (n <- n.) >= 0)      stopifnot(n == (n. <- as.integer(n)), (n <- n.) >= 0)
37      if((lx <- length(x)) == 1) x <- rep.int(x, n)      if((lx <- length(x)) == 1) x <- rep.int(x, n)
38      else if(lx != n) stop("length(x) must be 1 or n")      else if(lx != n) stop("length(x) must be 1 or n")
39      cls <-      stopifnot(is.character(shape), nchar(shape) == 1,
40          if(is.double(x)) "dsCMatrix"                any(shape == c("t","s","g"))) # triangular / symmetric / general
41          else if(is.logical(x)) "lsCMatrix"      kind <-
42            if(is.double(x)) "d"
43            else if(is.logical(x)) "l"
44          else { ## for now          else { ## for now
45              storage.mode(x) <- "double"              storage.mode(x) <- "double"
46              "dsCMatrix"              "d"
47          }          }
48      new(cls, Dim = c(n,n), x = x, uplo = uplo,      new(paste(kind, shape, "CMatrix", sep=''),
49            Dim = c(n,n), x = x, uplo = uplo,
50          i = if(n) 0:(n - 1L) else integer(0), p = 0:n)          i = if(n) 0:(n - 1L) else integer(0), p = 0:n)
51  }  }
52    
53    ## Pkg 'spdep' had (relatively slow) versions of this as_dsCMatrix_I()
54    .symDiagonal <- function(n, x = rep.int(1,n), uplo = "U")
55        .sparseDiagonal(n, x, uplo, shape = "s")
56    
57    ## instead of   diagU2N(as(Diagonal(n), "CsparseMatrix")), diag = "N" in any case:
58    .trDiagonal <- function(n, x = rep.int(1,n), uplo = "U")
59        .sparseDiagonal(n, x, uplo, shape = "t")
60    
61    
62  ### This is modified from a post of Bert Gunter to R-help on  1 Sep 2005.  ### This is modified from a post of Bert Gunter to R-help on  1 Sep 2005.
63  ### Bert's code built on a post by Andy Liaw who most probably was influenced  ### Bert's code built on a post by Andy Liaw who most probably was influenced
64  ### by earlier posts, notably one by Scott Chasalow on S-news, 16 Jan 2002  ### by earlier posts, notably one by Scott Chasalow on S-news, 16 Jan 2002
65  ### who posted his bdiag() function written in December 1995.  ### who posted his bdiag() function written in December 1995.
66    if(FALSE)##--- no longer used:
67  bdiag <- function(...) {  .bdiag <- function(lst) {
68      if(nargs() == 0) return(new("dgCMatrix"))      ### block-diagonal matrix [a dgTMatrix] from list of matrices
69      ## else :      stopifnot(is.list(lst), length(lst) >= 1)
70      mlist <- if (nargs() == 1) as.list(...) else list(...)      dims <- sapply(lst, dim, USE.NAMES=FALSE)
     dims <- sapply(mlist, dim)  
71      ## make sure we had all matrices:      ## make sure we had all matrices:
72      if(!(is.matrix(dims) && nrow(dims) == 2))      if(!(is.matrix(dims) && nrow(dims) == 2))
73          stop("some arguments are not matrices")          stop("some arguments are not matrices")
74      csdim <- rbind(rep.int(0L, 2),      csdim <- rbind(rep.int(0L, 2),
75                     apply(sapply(mlist, dim), 1, cumsum))                     apply(dims, 1, cumsum))
76      ret <- new("dgTMatrix", Dim = as.integer(csdim[nrow(csdim),]))      r <- new("dgTMatrix")
77        r@Dim <- as.integer(csdim[nrow(csdim),])
78      add1 <- matrix(1:0, 2,2)      add1 <- matrix(1:0, 2,2)
79      for(i in seq_along(mlist)) {      for(i in seq_along(lst)) {
80          indx <- apply(csdim[i:(i+1),] + add1, 2, function(n) n[1]:n[2])          indx <- apply(csdim[i:(i+1),] + add1, 2, function(n) n[1]:n[2])
81          if(is.null(dim(indx))) ## non-square matrix          if(is.null(dim(indx))) ## non-square matrix
82              ret[indx[[1]],indx[[2]]] <- mlist[[i]]              r[indx[[1]],indx[[2]]] <- lst[[i]]
83          else ## square matrix          else ## square matrix
84              ret[indx[,1],indx[,2]] <- mlist[[i]]              r[indx[,1], indx[,2]] <- lst[[i]]
85      }      }
86      ## slightly debatable if we really should return Csparse.. :      r
87      as(ret, "CsparseMatrix")  }
88    ## expand(<mer>) needed something like bdiag() for lower-triangular
89    ## (Tsparse) Matrices; hence Doug Bates provided a much more efficient
90    ##  implementation for those; now extended and generalized:
91    .bdiag <- function(lst) {
92        ## block-diagonal matrix [a dgTMatrix] from list of matrices
93        stopifnot(is.list(lst), (nl <- length(lst)) >= 1)
94    
95        Tlst <- lapply(lapply(lst, Matrix:::as_Csp2), # includes "diagU2N"
96                       as, "TsparseMatrix")
97        if(nl == 1) return(Tlst[[1]])
98        ## else
99        i_off <- c(0L, cumsum(sapply(Tlst, nrow)))
100        j_off <- c(0L, cumsum(sapply(Tlst, ncol)))
101    
102        clss <- sapply(Tlst, class)
103        knds <- substr(clss, 2, 2)
104        sym  <- knds == "s" # symmetric ones
105        tri  <- knds == "t" # triangular ones
106        use.n <- any(is.n <- substr(clss,1,1) == "n")
107        if(use.n && !(use.n <- all(is.n)))
108            Tlst[is.n] <- lapply(Tlst[is.n], as, "lMatrix")
109        if(all(sym)) { ## result should be *symmetric*
110            uplos <- sapply(Tlst, slot, "uplo") ## either "U" or "L"
111            tLU <- table(uplos)# of length 1 or 2 ..
112            if(length(tLU) == 1) { ## all "U" or all "L"
113                useU <- uplos[1] == "U"
114            } else { ## length(tLU) == 2, counting "L" and "U"
115                useU <- diff(tLU) >= 0
116                if(useU && (hasL <- tLU[1] > 0))
117                    Tlst[hasL] <- lapply(Tlst[hasL], t)
118                else if(!useU && (hasU <- tLU[2] > 0))
119                    Tlst[hasU] <- lapply(Tlst[hasU], t)
120            }
121            if(use.n) { ## return nsparseMatrix :
122                r <- new("nsTMatrix")
123            } else {
124                r <- new("dsTMatrix")
125                r@x <- unlist(lapply(Tlst, slot, "x"))
126            }
127            r@uplo <- if(useU) "U" else "L"
128        }
129        else if(all(tri) && { ULs <- sapply(Tlst, slot, "uplo")##  "U" or "L"
130                              all(ULs[1L] == ULs[-1L]) } ## all upper or all lower
131           ){ ## *triangular* result
132    
133            if(use.n) { ## return nsparseMatrix :
134                r <- new("ntTMatrix")
135            } else {
136                r <- new("dtTMatrix")
137                r@x <- unlist(lapply(Tlst, slot, "x"))
138            }
139            r@uplo <- ULs[1L]
140        }
141        else {
142            if(any(sym))
143                Tlst[sym] <- lapply(Tlst[sym], as, "generalMatrix")
144            if(use.n) { ## return nsparseMatrix :
145                r <- new("ngTMatrix")
146            } else {
147                r <- new("dgTMatrix")
148                r@x <- unlist(lapply(Tlst, slot, "x"))
149            }
150        }
151        r@Dim <- c(i_off[nl+1], j_off[nl + 1])
152        r@i <- unlist(lapply(1:nl, function(k) Tlst[[k]]@i + i_off[k]))
153        r@j <- unlist(lapply(1:nl, function(k) Tlst[[k]]@j + j_off[k]))
154        r
155    }
156    
157    bdiag <- function(...) {
158        if((nA <- nargs()) == 0) return(new("dgCMatrix"))
159        if(nA == 1 && !is.list(...))
160            return(as(..., "CsparseMatrix"))
161        alis <- if(nA == 1 && is.list(..1)) ..1 else list(...)
162        if(length(alis) == 1)
163            return(as(alis[[1]], "CsparseMatrix"))
164    
165        ## else : two or more arguments
166        as(.bdiag(alis), "CsparseMatrix")
167  }  }
168    
169    
# Line 325  Line 415 
415                                  j = "index", value = "replValue"),                                  j = "index", value = "replValue"),
416                   function(x,i,j, ..., value) replDiag(x, j=j, value=value))                   function(x,i,j, ..., value) replDiag(x, j=j, value=value))
417    
418    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",
419                                    value = "scarceMatrix"),
420                     function (x, i, j, ..., value)
421                     callGeneric(x=x, , j=j, value = as(value, "sparseVector")))
422    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "missing",
423                                    value = "scarceMatrix"),
424                     function (x, i, j, ..., value)
425                     callGeneric(x=x, i=i, , value = as(value, "sparseVector")))
426    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "index",
427                                    value = "scarceMatrix"),
428                     function (x, i, j, ..., value)
429                     callGeneric(x=x, i=i, j=j, value = as(value, "sparseVector")))
430    
431    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",
432                                    value = "sparseVector"),
433                     replDiag)
434    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "missing",
435                                    value = "sparseVector"),
436                     replDiag)
437    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "index",
438                                    value = "sparseVector"),
439                     replDiag)
440    
441    
442  setMethod("t", signature(x = "diagonalMatrix"),  setMethod("t", signature(x = "diagonalMatrix"),
443            function(x) { x@Dimnames <- x@Dimnames[2:1] ; x })            function(x) { x@Dimnames <- x@Dimnames[2:1] ; x })
# Line 619  Line 732 
732  ## ddi*: Arith: result numeric, potentially ddiMatrix  ## ddi*: Arith: result numeric, potentially ddiMatrix
733  setMethod("Arith", signature(e1 = "ddiMatrix", e2 = "numeric"),  setMethod("Arith", signature(e1 = "ddiMatrix", e2 = "numeric"),
734            function(e1,e2) {            function(e1,e2) {
735                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]; nsq <- n^2
736                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
737                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
738                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
# Line 637  Line 750 
750    
751  setMethod("Arith", signature(e1 = "numeric", e2 = "ddiMatrix"),  setMethod("Arith", signature(e1 = "numeric", e2 = "ddiMatrix"),
752            function(e1,e2) {            function(e1,e2) {
753                n <- e2@Dim[1]; nsq <- n*n                n <- e2@Dim[1]; nsq <- n^2
754                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
755                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
756                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
# Line 656  Line 769 
769  ## ldi* Arith --> result numeric, potentially ddiMatrix  ## ldi* Arith --> result numeric, potentially ddiMatrix
770  setMethod("Arith", signature(e1 = "ldiMatrix", e2 = "numeric"),  setMethod("Arith", signature(e1 = "ldiMatrix", e2 = "numeric"),
771            function(e1,e2) {            function(e1,e2) {
772                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]; nsq <- n^2
773                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
774                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
775                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
# Line 675  Line 788 
788    
789  setMethod("Arith", signature(e1 = "numeric", e2 = "ldiMatrix"),  setMethod("Arith", signature(e1 = "numeric", e2 = "ldiMatrix"),
790            function(e1,e2) {            function(e1,e2) {
791                n <- e2@Dim[1]; nsq <- n*n                n <- e2@Dim[1]; nsq <- n^2
792                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
793                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
794                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
# Line 695  Line 808 
808  ## ddi*: for "Ops" without Arith --> result logical, potentially ldi  ## ddi*: for "Ops" without Arith --> result logical, potentially ldi
809  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "numeric"),  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "numeric"),
810            function(e1,e2) {            function(e1,e2) {
811                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]; nsq <- n^2
812                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
813                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
814                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
# Line 714  Line 827 
827    
828  setMethod("Ops", signature(e1 = "numeric", e2 = "ddiMatrix"),  setMethod("Ops", signature(e1 = "numeric", e2 = "ddiMatrix"),
829            function(e1,e2) {            function(e1,e2) {
830                n <- e2@Dim[1]; nsq <- n*n                n <- e2@Dim[1]; nsq <- n^2
831                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
832                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
833                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
# Line 734  Line 847 
847  ## ldi*: for "Ops" without Arith --> result logical, potentially ldi  ## ldi*: for "Ops" without Arith --> result logical, potentially ldi
848  setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "numeric"),  setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "numeric"),
849            function(e1,e2) {            function(e1,e2) {
850                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]; nsq <- n^2
851                f0 <- callGeneric(FALSE, e2)                f0 <- callGeneric(FALSE, e2)
852                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
853                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
# Line 752  Line 865 
865    
866  setMethod("Ops", signature(e1 = "numeric", e2 = "ldiMatrix"),  setMethod("Ops", signature(e1 = "numeric", e2 = "ldiMatrix"),
867            function(e1,e2) {            function(e1,e2) {
868                n <- e2@Dim[1]; nsq <- n*n                n <- e2@Dim[1]; nsq <- n^2
869                f0 <- callGeneric(e1, FALSE)                f0 <- callGeneric(e1, FALSE)
870                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
871                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
# Line 851  Line 964 
964            function(object) {            function(object) {
965                d <- dim(object)                d <- dim(object)
966                cl <- class(object)                cl <- class(object)
967                cat(sprintf('%d x %d diagonal matrix of class "%s"\n',                cat(sprintf('%d x %d diagonal matrix of class "%s"',
968                            d[1], d[2], cl))                            d[1], d[2], cl))
969                  if(d[1] < 50) {
970                      cat("\n")
971                prDiag(object)                prDiag(object)
972                  } else {
973                      cat(", with diagonal entries\n")
974                      show(diag(object))
975                      invisible(object)
976                  }
977            })            })

Legend:
Removed from v.2197  
changed lines
  Added in v.2226

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