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 1635, Tue Oct 17 21:21:44 2006 UTC revision 1654, Fri Oct 27 16:58:15 2006 UTC
# Line 13  Line 13 
13          n <- as.integer(n)          n <- as.integer(n)
14      }      }
15    
16      if(missing(x)) # unit diagonal matrix      if(missing(x)) ## unit diagonal matrix
17          new("ddiMatrix", Dim = c(n,n), diag = "U")          new("ddiMatrix", Dim = c(n,n), diag = "U")
18      else {      else {
19          stopifnot(length(x) == n)          stopifnot(length(x) == n)
# Line 47  Line 47 
47                     apply(sapply(mlist, dim), 1, cumsum))                     apply(sapply(mlist, dim), 1, cumsum))
48      ret <- new("dgTMatrix", Dim = as.integer(csdim[nrow(csdim),]))      ret <- new("dgTMatrix", Dim = as.integer(csdim[nrow(csdim),]))
49      add1 <- matrix(1:0, 2,2)      add1 <- matrix(1:0, 2,2)
50      for(i in seq(along = mlist)) {      for(i in seq_along(mlist)) {
51          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])
52          if(is.null(dim(indx))) ## non-square matrix          if(is.null(dim(indx))) ## non-square matrix
53              ret[indx[[1]],indx[[2]]] <- mlist[[i]]              ret[indx[[1]],indx[[2]]] <- mlist[[i]]
# Line 58  Line 58 
58      as(ret, "CsparseMatrix")      as(ret, "CsparseMatrix")
59  }  }
60    
61  setAs("diagonalMatrix", "triangularMatrix",  diag2T <- function(from) {
62        function(from) {      i <- if(from@diag == "U") integer(0) else seq_len(from@Dim[1]) - 1:1
           n <- from@Dim[1]  
           i <- seq(length = n)  
           x <- from@x  
63            new(paste(.M.kind(from), "tTMatrix", sep=''),            new(paste(.M.kind(from), "tTMatrix", sep=''),
64                diag = from@diag, Dim = from@Dim, Dimnames = from@Dimnames,                diag = from@diag, Dim = from@Dim, Dimnames = from@Dimnames,
65                x = x, i = i, j = i)          x = from@x, # <- ok for diag = "U" and "N" (!)
66            })          i = i, j = i)
67    }
68    
69    setAs("diagonalMatrix", "triangularMatrix", diag2T)
70    setAs("diagonalMatrix", "sparseMatrix", diag2T)
71    ## is better than this:
72    ## setAs("diagonalMatrix", "sparseMatrix",
73    ##       function(from)
74    ##        as(from, if(is(from, "dMatrix")) "dgCMatrix" else "lgCMatrix"))
75    setAs("diagonalMatrix", "CsparseMatrix",
76          function(from) as(diag2T(from), "CsparseMatrix"))
77    
78  setAs("diagonalMatrix", "matrix",  setAs("diagonalMatrix", "matrix",
79        function(from) {        function(from) {
# Line 76  Line 83 
83                 nrow = n, ncol = n)                 nrow = n, ncol = n)
84        })        })
85    
86  setAs("diagonalMatrix", "generalMatrix",  setAs("diagonalMatrix", "generalMatrix", # prefer sparse:
87        function(from) as(from, paste(.M.kind(from), "geMatrix", sep='')))        function(from) as(from, paste(.M.kind(from), "gCMatrix", sep='')))
88    
89    ## given the above, the following  4  coercions should be all unneeded;
90    ## we prefer triangular to general:
91  setAs("ddiMatrix", "dgTMatrix",  setAs("ddiMatrix", "dgTMatrix",
92        function(from) {        function(from) {
93              .Deprecated("as(, \"sparseMatrix\")")
94            n <- from@Dim[1]            n <- from@Dim[1]
95            i <- seq(length = n) - 1:1            i <- seq_len(n) - 1:1
96            new("dgTMatrix", i = i, j = i,            new("dgTMatrix", i = i, j = i,
97                x = if(from@diag == "U") rep(1,n) else from@x,                x = if(from@diag == "U") rep(1,n) else from@x,
98                Dim = c(n,n), Dimnames = from@Dimnames) })                Dim = c(n,n), Dimnames = from@Dimnames) })
# Line 92  Line 102 
102    
103  setAs("ldiMatrix", "lgTMatrix",  setAs("ldiMatrix", "lgTMatrix",
104        function(from) {        function(from) {
105              .Deprecated("as(, \"sparseMatrix\")")
106            n <- from@Dim[1]            n <- from@Dim[1]
107            if(from@diag == "U") { # unit-diagonal            if(from@diag == "U") { # unit-diagonal
108                x <- rep.int(TRUE, n)                x <- rep.int(TRUE, n)
109                i <- seq(length = n)                i <- seq_len(n) - 1:1
110            } else { # "normal"            } else { # "normal"
111                nz <- nz.NA(from@x, na. = TRUE)                nz <- nz.NA(from@x, na. = TRUE)
112                x <- from@x[nz]                x <- from@x[nz]
# Line 107  Line 118 
118  setAs("ldiMatrix", "lgCMatrix",  setAs("ldiMatrix", "lgCMatrix",
119        function(from) as(as(from, "lgTMatrix"), "lgCMatrix"))        function(from) as(as(from, "lgTMatrix"), "lgCMatrix"))
120    
 setAs("diagonalMatrix", "sparseMatrix",  
       function(from)  
           as(from, if(is(from, "dMatrix")) "dgCMatrix" else "lgCMatrix"))  
121    
122  if(FALSE) # now have faster  "ddense" -> "dge"  if(FALSE) # now have faster  "ddense" -> "dge"
123  setAs("ddiMatrix", "dgeMatrix",  setAs("ddiMatrix", "dgeMatrix",
# Line 175  Line 183 
183  setMethod("isSymmetric", signature(object = "diagonalMatrix"),  setMethod("isSymmetric", signature(object = "diagonalMatrix"),
184            function(object) TRUE)            function(object) TRUE)
185    
186    setMethod("chol", signature(x = "ddiMatrix"),# pivot = "ANY"
187              function(x, pivot) {
188                  if(any(x@x < 0)) stop("chol() is undefined for diagonal matrix with negative entries")
189                  x@x <- sqrt(x@x)
190                  x
191              })
192    ## chol(L) is L for logical diagonal:
193    setMethod("chol", signature(x = "ldiMatrix"), function(x, pivot) x)
194    
195    
196  setMethod("diag", signature(x = "diagonalMatrix"),  setMethod("diag", signature(x = "diagonalMatrix"),
197            function(x = 1, nrow, ncol = n) {            function(x = 1, nrow, ncol = n) {
198               if(x@diag == "U")               if(x@diag == "U")
# Line 193  Line 211 
211  })  })
212    
213  ## Basic Matrix Multiplication {many more to add}  ## Basic Matrix Multiplication {many more to add}
214    ##       ---------------------
215  ## FIXME: extend this for 'ldi', i.e. do "diagonalMatrix"  ## Note that "ldi" logical are treated as numeric
216  diagdiagprod <- function(x, y) {  diagdiagprod <- function(x, y) {
217      if(any(dim(x) != dim(y))) stop("non-matching dimensions")      if(any(dim(x) != dim(y))) stop("non-matching dimensions")
218      if(x@diag != "U") {      if(x@diag != "U") {
219          if(y@diag != "U") x@x <- x@x * y@x          if(y@diag != "U") {
220                nx <- x@x * y@x
221                if(is.numeric(nx) && !is.numeric(x@x))
222                    x <- as(x, "dMatrix")
223                x@x <- as.numeric(nx)
224            }
225          return(x)          return(x)
226      } else ## x is unit diagonal      } else ## x is unit diagonal
227      return(y)      return(y)
228  }  }
229    
230  setMethod("%*%", signature(x = "ddiMatrix", y = "ddiMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "diagonalMatrix"),
231            diagdiagprod, valueClass = "ddiMatrix")            diagdiagprod, valueClass = "ddiMatrix")
232    
233  formals(diagdiagprod) <- alist(x=, y=NULL)  formals(diagdiagprod) <- alist(x=, y=x)
234  setMethod("crossprod", signature(x = "ddiMatrix", y = "ddiMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "diagonalMatrix"),
235              diagdiagprod, valueClass = "ddiMatrix")
236    setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "diagonalMatrix"),
237            diagdiagprod, valueClass = "ddiMatrix")            diagdiagprod, valueClass = "ddiMatrix")
238  setMethod("tcrossprod", signature(x = "ddiMatrix", y = "ddiMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "missing"),
239              diagdiagprod, valueClass = "ddiMatrix")
240    setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "missing"),
241            diagdiagprod, valueClass = "ddiMatrix")            diagdiagprod, valueClass = "ddiMatrix")
242    
243    
# Line 227  Line 254 
254  formals(diagmatprod) <- alist(x=, y=NULL)  formals(diagmatprod) <- alist(x=, y=NULL)
255  setMethod("crossprod", signature(x = "diagonalMatrix", y = "matrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "matrix"),
256            diagmatprod)            diagmatprod)
257    setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "matrix"),
258              diagmatprod)
259    
260  diagdgeprod <- function(x, y) {  diagdgeprod <- function(x, y) {
261      dx <- dim(x)      dx <- dim(x)
# Line 278  Line 307 
307    
308  ## These are cheap implementations via coercion  ## These are cheap implementations via coercion
309    
310    setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "sparseMatrix"),
311              function(e1,e2) callGeneric(as(e1, "sparseMatrix"), e2))
312    setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "diagonalMatrix"),
313              function(e1,e2) callGeneric(e1, as(e2, "sparseMatrix")))
314    
315  ## FIXME?: In theory, this can be done *FASTER*, in some cases, via tapply1()  ## FIXME?: In theory, this can be done *FASTER*, in some cases, via tapply1()
316    
317  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),

Legend:
Removed from v.1635  
changed lines
  Added in v.1654

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