SCM Repository

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

Diff of /pkg/R/diagMatrix.R

revision 1653, Tue Oct 24 21:45:22 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),]))
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.1653 changed lines Added in v.1654