SCM Repository

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

Diff of /pkg/R/diagMatrix.R

revision 1635, Tue Oct 17 21:21:44 2006 UTC revision 1655, Mon Oct 30 17:16:27 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    .diag.x <- function(m) {
90        if(m@diag == "U")
91            rep.int(if(is.numeric(m@x)) 1. else TRUE,
92                    m@Dim[1])
93        else m@x
94    }
95
96    .diag.2N <- function(m) {
97        if(m@diag == "U") m@diag <- "N"
98        m
99    }
100
101    ## given the above, the following  4  coercions should be all unneeded;
102    ## we prefer triangular to general:
103  setAs("ddiMatrix", "dgTMatrix",  setAs("ddiMatrix", "dgTMatrix",
104        function(from) {        function(from) {
105              .Deprecated("as(, \"sparseMatrix\")")
106            n <- from@Dim[1]            n <- from@Dim[1]
107            i <- seq(length = n) - 1:1            i <- seq_len(n) - 1:1
108            new("dgTMatrix", i = i, j = i,            new("dgTMatrix", i = i, j = i, x = .diag.x(from),
x = if(from@diag == "U") rep(1,n) else from@x,
109                Dim = c(n,n), Dimnames = from@Dimnames) })                Dim = c(n,n), Dimnames = from@Dimnames) })
110
111  setAs("ddiMatrix", "dgCMatrix",  setAs("ddiMatrix", "dgCMatrix",
112        function(from) as(as(from, "dgTMatrix"), "dgCMatrix"))        function(from) as(as(from, "sparseMatrix"), "dgCMatrix"))
113
114  setAs("ldiMatrix", "lgTMatrix",  setAs("ldiMatrix", "lgTMatrix",
115        function(from) {        function(from) {
116              .Deprecated("as(, \"sparseMatrix\")")
117            n <- from@Dim[1]            n <- from@Dim[1]
118            if(from@diag == "U") { # unit-diagonal            if(from@diag == "U") { # unit-diagonal
119                x <- rep.int(TRUE, n)                x <- rep.int(TRUE, n)
120                i <- seq(length = n)                i <- seq_len(n) - 1:1
121            } else { # "normal"            } else { # "normal"
122                nz <- nz.NA(from@x, na. = TRUE)                nz <- nz.NA(from@x, na. = TRUE)
123                x <- from@x[nz]                x <- from@x[nz]
# Line 107  Line 129
129  setAs("ldiMatrix", "lgCMatrix",  setAs("ldiMatrix", "lgCMatrix",
130        function(from) as(as(from, "lgTMatrix"), "lgCMatrix"))        function(from) as(as(from, "lgTMatrix"), "lgCMatrix"))
131
setAs("diagonalMatrix", "sparseMatrix",
function(from)
as(from, if(is(from, "dMatrix")) "dgCMatrix" else "lgCMatrix"))
132
133  if(FALSE) # now have faster  "ddense" -> "dge"  if(FALSE) # now have faster  "ddense" -> "dge"
134  setAs("ddiMatrix", "dgeMatrix",  setAs("ddiMatrix", "dgeMatrix",
# Line 175  Line 194
194  setMethod("isSymmetric", signature(object = "diagonalMatrix"),  setMethod("isSymmetric", signature(object = "diagonalMatrix"),
195            function(object) TRUE)            function(object) TRUE)
196
197    setMethod("chol", signature(x = "ddiMatrix"),# pivot = "ANY"
198              function(x, pivot) {
199                  if(any(x@x < 0)) stop("chol() is undefined for diagonal matrix with negative entries")
200                  x@x <- sqrt(x@x)
201                  x
202              })
203    ## chol(L) is L for logical diagonal:
204    setMethod("chol", signature(x = "ldiMatrix"), function(x, pivot) x)
205
206
207  setMethod("diag", signature(x = "diagonalMatrix"),  setMethod("diag", signature(x = "diagonalMatrix"),
208            function(x = 1, nrow, ncol = n) {            function(x = 1, nrow, ncol = n) {
209               if(x@diag == "U")               if(x@diag == "U")
# Line 193  Line 222
222  })  })
223
224  ## Basic Matrix Multiplication {many more to add}  ## Basic Matrix Multiplication {many more to add}
225    ##       ---------------------
226  ## FIXME: extend this for 'ldi', i.e. do "diagonalMatrix"  ## Note that "ldi" logical are treated as numeric
227  diagdiagprod <- function(x, y) {  diagdiagprod <- function(x, y) {
228      if(any(dim(x) != dim(y))) stop("non-matching dimensions")      if(any(dim(x) != dim(y))) stop("non-matching dimensions")
229      if(x@diag != "U") {      if(x@diag != "U") {
230          if(y@diag != "U") x@x <- x@x * y@x          if(y@diag != "U") {
231                nx <- x@x * y@x
232                if(is.numeric(nx) && !is.numeric(x@x))
233                    x <- as(x, "dMatrix")
234                x@x <- as.numeric(nx)
235            }
236          return(x)          return(x)
237      } else ## x is unit diagonal      } else ## x is unit diagonal
238      return(y)      return(y)
239  }  }
240
241  setMethod("%*%", signature(x = "ddiMatrix", y = "ddiMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "diagonalMatrix"),
242            diagdiagprod, valueClass = "ddiMatrix")            diagdiagprod, valueClass = "ddiMatrix")
243
244  formals(diagdiagprod) <- alist(x=, y=NULL)  formals(diagdiagprod) <- alist(x=, y=x)
245  setMethod("crossprod", signature(x = "ddiMatrix", y = "ddiMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "diagonalMatrix"),
246            diagdiagprod, valueClass = "ddiMatrix")            diagdiagprod, valueClass = "ddiMatrix")
247  setMethod("tcrossprod", signature(x = "ddiMatrix", y = "ddiMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "diagonalMatrix"),
248              diagdiagprod, valueClass = "ddiMatrix")
249    setMethod("crossprod", signature(x = "diagonalMatrix", y = "missing"),
250              diagdiagprod, valueClass = "ddiMatrix")
251    setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "missing"),
252            diagdiagprod, valueClass = "ddiMatrix")            diagdiagprod, valueClass = "ddiMatrix")
253
254
# Line 227  Line 265
265  formals(diagmatprod) <- alist(x=, y=NULL)  formals(diagmatprod) <- alist(x=, y=NULL)
266  setMethod("crossprod", signature(x = "diagonalMatrix", y = "matrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "matrix"),
267            diagmatprod)            diagmatprod)
268    setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "matrix"),
269              diagmatprod)
270
271  diagdgeprod <- function(x, y) {  diagdgeprod <- function(x, y) {
272      dx <- dim(x)      dx <- dim(x)
# Line 276  Line 316
316
317  ### ---------------- diagonal  o   sparse  -----------------------------  ### ---------------- diagonal  o   sparse  -----------------------------
318
319
320    ## Use function for several signatures, in order to evade
321    ## ambiguous dispatch for "ddi", since there's also Arith(ddense., ddense.)
322    diagOdiag <- function(e1,e2) { # result should also be diagonal
323        r <- callGeneric(.diag.x(e1), .diag.x(e2)) # error if not "compatible"
324        if(is.numeric(r)) {
325            if(is.numeric(e2@x)) {
326                e2@x <- r; return(.diag.2N(e2)) }
327            if(!is.numeric(e1@x))
328                ## e.g. e1, e2 are logical;
329                e1 <- as(e1, "dMatrix")
330        }
331        else if(is.logical(r))
332            e1 <- as(e1, "lMatrix")
333        else stop("intermediate 'r' is of type", typeof(r))
334        e1@x <- r
335        .diag.2N(e1)
336    }
337
338    setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "diagonalMatrix"),
339              diagOdiag)
340    ## These two are just for method disambiguation:
341    setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "diagonalMatrix"),
342              diagOdiag)
343    setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "ddiMatrix"),
344              diagOdiag)
345
346    ## For almost everything else, diag* shall be treated "as sparse" :
347  ## These are cheap implementations via coercion  ## These are cheap implementations via coercion
348
349    ## for disambiguation
350    setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "sparseMatrix"),
351              function(e1,e2) callGeneric(as(e1, "sparseMatrix"), e2))
352    setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "diagonalMatrix"),
353              function(e1,e2) callGeneric(e1, as(e2, "sparseMatrix")))
354    ## in general:
355    setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "ANY"),
356              function(e1,e2) callGeneric(as(e1,"sparseMatrix"), e2))
357    setMethod("Ops", signature(e1 = "ANY", e2 = "diagonalMatrix"),
358              function(e1,e2) callGeneric(e1, as(e2,"sparseMatrix")))
359
360
361
362  ## 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()
363
364  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),

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