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 2052, Wed Aug 15 13:33:19 2007 UTC revision 2185, Sat Apr 26 20:33:16 2008 UTC
# Line 16  Line 16 
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)          lx <- length(x)
20            stopifnot(lx == 1 || lx == n) # but keep 'x' short for now
21          if(is.logical(x))          if(is.logical(x))
22              cl <- "ldiMatrix"              cl <- "ldiMatrix"
23          else if(is.numeric(x)) {          else if(is.numeric(x)) {
# Line 26  Line 27 
27          else if(is.complex(x)) {          else if(is.complex(x)) {
28              cl <- "zdiMatrix"  # will not yet work              cl <- "zdiMatrix"  # will not yet work
29          } else stop("'x' has invalid data type")          } else stop("'x' has invalid data type")
30          new(cl, Dim = c(n,n), diag = "N", x = x)          new(cl, Dim = c(n,n), diag = "N",
31                x = if(lx == 1) rep.int(x,n) else x)
32      }      }
33  }  }
34    
35    ## Pkg 'spdep' had (relatively slow) versions of this as_dsCMatrix_I()
36    .symDiagonal <- function(n, x = rep.int(1,n), uplo = "U") {
37        stopifnot(n == (n. <- as.integer(n)), (n <- n.) >= 0)
38        if((lx <- length(x)) == 1) x <- rep.int(x, n)
39        else if(lx != n) stop("length(x) must be 1 or n")
40        cls <-
41            if(is.double(x)) "dsCMatrix"
42            else if(is.logical(x)) "lsCMatrix"
43            else { ## for now
44                storage.mode(x) <- "double"
45                "dsCMatrix"
46            }
47        new(cls, Dim = c(n,n), x = x, uplo = uplo,
48            i = if(n) 0:(n - 1L) else integer(0), p = 0:n)
49    }
50    
51  ### 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.
52  ### 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
53  ### 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
# Line 58  Line 76 
76      as(ret, "CsparseMatrix")      as(ret, "CsparseMatrix")
77  }  }
78    
79  diag2tT <- function(from) {  
80    .diag2tT <- function(from, uplo = "U", kind = .M.kind(from)) {
81        ## to triangular Tsparse
82      i <- if(from@diag == "U") integer(0) else seq_len(from@Dim[1]) - 1L      i <- if(from@diag == "U") integer(0) else seq_len(from@Dim[1]) - 1L
83      new(paste(.M.kind(from), "tTMatrix", sep=''),      new(paste(kind, "tTMatrix", sep=''),
84          diag = from@diag, Dim = from@Dim, Dimnames = from@Dimnames,          diag = from@diag, Dim = from@Dim, Dimnames = from@Dimnames,
85            uplo = uplo,
86          x = from@x, # <- ok for diag = "U" and "N" (!)          x = from@x, # <- ok for diag = "U" and "N" (!)
87          i = i, j = i)          i = i, j = i)
88  }  }
89    
90  diag2sT <- function(from) { # to symmetric Tsparse  .diag2sT <- function(from, uplo = "U", kind = .M.kind(from)) {
91      i <- if(from@diag == "U") integer(0) else seq_len(from@Dim[1]) - 1L      ## to symmetric Tsparse
92      new(paste(.M.kind(from), "sTMatrix", sep=''),      n <- from@Dim[1]
93        i <- seq_len(n) - 1L
94        new(paste(kind, "sTMatrix", sep=''),
95          Dim = from@Dim, Dimnames = from@Dimnames,          Dim = from@Dim, Dimnames = from@Dimnames,
96          x = from@x, i = i, j = i)          i = i, j = i, uplo = uplo,
97            x = if(from@diag == "N") from@x else ## "U"-diag
98            rep.int(switch(kind,
99                           "d" = 1.,
100                           "l" =,
101                           "n" = TRUE,
102                           ## otherwise
103                           stop("'", kind,"' kind not yet implemented")), n))
104    }
105    
106    ## diagonal -> triangular,  upper / lower depending on "partner":
107    diag2tT.u <- function(d, x, kind = .M.kind(d))
108        .diag2tT(d, uplo = if(is(x,"triangularMatrix")) x@uplo else "U", kind)
109    
110    ## diagonal -> sparse {triangular OR symmetric} (upper / lower) depending on "partner":
111    diag2Tsmart <- function(d, x, kind = .M.kind(d)) {
112        clx <- getClassDef(class(x))
113        if(extends(clx, "symmetricMatrix"))
114            .diag2sT(d, uplo = x@uplo, kind)
115        else
116            .diag2tT(d, uplo = if(extends(clx,"triangularMatrix")) x@uplo else "U", kind)
117  }  }
118    
119  setAs("diagonalMatrix", "triangularMatrix", diag2tT)  
120  setAs("diagonalMatrix", "sparseMatrix", diag2tT)  ## In order to evade method dispatch ambiguity warnings,
121    ## and because we can save a .M.kind() call, we use this explicit
122    ## "hack"  instead of signature  x = "diagonalMatrix" :
123    ##
124    ## ddi*:
125    diag2tT <- function(from) .diag2tT(from, "U", "d")
126    setAs("ddiMatrix", "triangularMatrix", diag2tT)
127    setAs("ddiMatrix", "sparseMatrix", diag2tT)
128  ## needed too (otherwise <dense> -> Tsparse is taken):  ## needed too (otherwise <dense> -> Tsparse is taken):
129  setAs("diagonalMatrix", "TsparseMatrix", diag2tT)  setAs("ddiMatrix", "TsparseMatrix", diag2tT)
130  ## is better than this:  setAs("ddiMatrix", "CsparseMatrix",
131  ## setAs("diagonalMatrix", "sparseMatrix",        function(from) as(.diag2tT(from, "U", "d"), "CsparseMatrix"))
132  ##       function(from)  setAs("ddiMatrix", "symmetricMatrix",
133  ##        as(from, if(is(from, "dMatrix")) "dgCMatrix" else "lgCMatrix"))        function(from) .diag2sT(from, "U", "d"))
134  setAs("diagonalMatrix", "CsparseMatrix",  ##
135        function(from) as(diag2tT(from), "CsparseMatrix"))  ## ldi*:
136    diag2tT <- function(from) .diag2tT(from, "U", "l")
137    setAs("ldiMatrix", "triangularMatrix", diag2tT)
138    setAs("ldiMatrix", "sparseMatrix", diag2tT)
139    ## needed too (otherwise <dense> -> Tsparse is taken):
140    setAs("ldiMatrix", "TsparseMatrix", diag2tT)
141    setAs("ldiMatrix", "CsparseMatrix",
142          function(from) as(.diag2tT(from, "U", "l"), "CsparseMatrix"))
143    setAs("ldiMatrix", "symmetricMatrix",
144          function(from) .diag2sT(from, "U", "l"))
145    
 setAs("diagonalMatrix", "symmetricMatrix", diag2sT)  
146    
147  setAs("diagonalMatrix", "matrix",  setAs("diagonalMatrix", "nMatrix",
148        function(from) {        function(from) {
149            n <- from@Dim[1]            n <- from@Dim[1]
150            diag(x = if(from@diag == "U") { if(is.numeric(from@x)) 1. else TRUE            i <- if(from@diag == "U") integer(0) else which(isN0(from@x)) - 1L
151                                       } else from@x,            new("ntTMatrix", i = i, j = i, diag = from@diag,
152                 nrow = n, ncol = n)                Dim = from@Dim, Dimnames = from@Dimnames)
153        })        })
154    
 setAs("diagonalMatrix", "generalMatrix", # prefer sparse:  
       function(from) as(from, paste(.M.kind(from), "gCMatrix", sep='')))  
   
 .diag.x <- function(m) {  
     if(m@diag == "U")  
         rep.int(if(is.numeric(m@x)) 1. else TRUE,  
                 m@Dim[1])  
     else m@x  
 }  
155    
156  .diag.2N <- function(m) {  ## Cheap fast substitute for diag() which *does* preserve the mode of x :
157      if(m@diag == "U") m@diag <- "N"  mkDiag <- function(x, n) {
158      m      y <- matrix(as0(mod=mode(x)), n,n)
159        if (n > 0) y[1L + 0:(n - 1L) * (n + 1L)] <- x
160        y
161  }  }
162    
163  ## given the above, the following  4  coercions should be all unneeded;  setAs("diagonalMatrix", "matrix",
 ## we prefer triangular to general:  
 setAs("ddiMatrix", "dgTMatrix",  
164        function(from) {        function(from) {
165            .Deprecated("as(, \"sparseMatrix\")")            ## want "ldiMatrix" -> <logical> "matrix" :
166            n <- from@Dim[1]            mkDiag(if(from@diag == "U") as1(from@x) else from@x,
167            i <- seq_len(n) - 1L                   n = from@Dim[1])
168            new("dgTMatrix", i = i, j = i, x = .diag.x(from),        })
               Dim = c(n,n), Dimnames = from@Dimnames) })  
169    
170  setAs("ddiMatrix", "dgCMatrix",  setMethod("as.vector", signature(x = "diagonalMatrix", mode="missing"),
171        function(from) as(as(from, "sparseMatrix"), "dgCMatrix"))            function(x, mode) {
172                  n <- x@Dim[1]
173                  mod.x <- mode(x@x)
174                  r <- vector(mod.x, length = n^2)
175                  if(n)
176                      r[1 + 0:(n - 1) * (n + 1)] <-
177                          if(x@diag == "U") as1(mod=mod.x) else x@x
178                  r
179              })
180    
181  setAs("ldiMatrix", "lgTMatrix",  setAs("diagonalMatrix", "generalMatrix", # prefer sparse:
182        function(from) {        function(from) as(as(from, "CsparseMatrix"), "generalMatrix"))
           .Deprecated("as(, \"sparseMatrix\")")  
           n <- from@Dim[1]  
           if(from@diag == "U") { # unit-diagonal  
               x <- rep.int(TRUE, n)  
               i <- seq_len(n) - 1L  
           } else { # "normal"  
               nz <- nz.NA(from@x, na. = TRUE)  
               x <- from@x[nz]  
               i <- which(nz) - 1L  
           }  
           new("lgTMatrix", i = i, j = i, x = x,  
               Dim = c(n,n), Dimnames = from@Dimnames) })  
183    
184  setAs("ldiMatrix", "lgCMatrix",  .diag.x <- function(m) if(m@diag == "U") rep.int(as1(m@x), m@Dim[1]) else m@x
       function(from) as(as(from, "lgTMatrix"), "lgCMatrix"))  
185    
186    .diag.2N <- function(m) {
187        if(m@diag == "U") m@diag <- "N"
188        m
189    }
190    
 if(FALSE) # now have faster  "ddense" -> "dge"  
191  setAs("ddiMatrix", "dgeMatrix",  setAs("ddiMatrix", "dgeMatrix",
192        function(from) as(as(from, "matrix"), "dgeMatrix"))        function(from) .Call(dup_mMatrix_as_dgeMatrix, from))
193    setAs("ddiMatrix", "ddenseMatrix",
194          function(from) as(as(from, "triangularMatrix"),"denseMatrix"))
195    setAs("ldiMatrix", "ldenseMatrix",
196          function(from) as(as(from, "triangularMatrix"),"denseMatrix"))
197    
198    
199  setAs("matrix", "diagonalMatrix",  setAs("matrix", "diagonalMatrix",
200        function(from) {        function(from) {
# Line 188  Line 239 
239  setMethod("diag", signature(x = "diagonalMatrix"),  setMethod("diag", signature(x = "diagonalMatrix"),
240            function(x = 1, nrow, ncol) .diag.x(x))            function(x = 1, nrow, ncol) .diag.x(x))
241    
242    subDiag <- function(x, i, j, ..., drop) {
 subDiag <- function(x, i, j, drop) {  
243      x <- as(x, "sparseMatrix")      x <- as(x, "sparseMatrix")
244      x <- if(missing(i))      x <- if(missing(i))
245          x[, j, drop=drop]          x[, j, drop=drop]
# Line 197  Line 247 
247          x[i, , drop=drop]          x[i, , drop=drop]
248      else      else
249          x[i,j, drop=drop]          x[i,j, drop=drop]
250      if(isDiagonal(x)) as(x, "diagonalMatrix") else x      if(isS4(x) && isDiagonal(x)) as(x, "diagonalMatrix") else x
251  }  }
252    
253  setMethod("[", signature(x = "diagonalMatrix", i = "index",  setMethod("[", signature(x = "diagonalMatrix", i = "index",
254                           j = "index", drop = "logical"), subDiag)                           j = "index", drop = "logical"), subDiag)
255  setMethod("[", signature(x = "diagonalMatrix", i = "index",  setMethod("[", signature(x = "diagonalMatrix", i = "index",
256                          j = "missing", drop = "logical"),                          j = "missing", drop = "logical"),
257            function(x, i, drop) subDiag(x, i=i, drop=drop))            function(x, i, j, ..., drop) subDiag(x, i=i, drop=drop))
258  setMethod("[", signature(x = "diagonalMatrix", i = "missing",  setMethod("[", signature(x = "diagonalMatrix", i = "missing",
259                           j = "index", drop = "logical"),                           j = "index", drop = "logical"),
260            function(x, j, drop) subDiag(x, j=j, drop=drop))            function(x, i, j, ..., drop) subDiag(x, j=j, drop=drop))
261    
262  ## When you assign to a diagonalMatrix, the result should be  ## When you assign to a diagonalMatrix, the result should be
263  ## diagonal or sparse ---  ## diagonal or sparse ---
264  ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch  ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch
265  replDiag <- function(x, i, j, value) {  ## Only(?) current bug:  x[i] <- value  is wrong when  i is *vector*
266    replDiag <- function(x, i, j, ..., value) {
267      x <- as(x, "sparseMatrix")      x <- as(x, "sparseMatrix")
268      if(missing(i))      if(missing(i))
269          x[, j] <- value          x[, j] <- value
270      else if(missing(j))      else if(missing(j)) { ##  x[i , ] <- v  *OR*   x[i] <- v
271            na <- nargs()
272    ##         message("diagnosing replDiag() -- nargs()= ", na)
273            if(na == 4)
274          x[i, ] <- value          x[i, ] <- value
275      else          else if(na == 3)
276                x[i] <- value
277            else stop("Internal bug: nargs()=",na,"; please report")
278        } else
279          x[i,j] <- value          x[i,j] <- value
280      if(isDiagonal(x)) as(x, "diagonalMatrix") else x      if(isDiagonal(x)) as(x, "diagonalMatrix") else x
281  }  }
282    
283  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index",
284                                  j = "index", value = "replValue"), replDiag)                                  j = "index", value = "replValue"), replDiag)
285    
286  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index",
287                                  j = "missing", value = "replValue"),                                  j = "missing", value = "replValue"),
288                   function(x, i, value) replDiag(x, i=i, value=value))                   function(x,i,j, ..., value) {
289                         ## message("before replDiag() -- nargs()= ", nargs())
290                         if(nargs() == 3)
291                             replDiag(x, i=i, value=value)
292                         else ## nargs() == 4 :
293                             replDiag(x, i=i, , value=value)
294                     })
295    
296    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "matrix", # 2-col.matrix
297                                    j = "missing", value = "replValue"),
298                     function(x,i,j, ..., value) {
299                         if(ncol(i) == 2) {
300                             if(all((ii <- i[,1]) == i[,2])) { # replace in diagonal only
301                                 x@x[ii] <- value
302                                 x
303                             } else { ## no longer diagonal, but remain sparse:
304                                 x <- as(x, "sparseMatrix")
305                                 x[i] <- value
306                                 x
307                             }
308                         }
309                         else { # behave as "base R": use as if vector
310                             x <- as(x, "matrix")
311                             x[i] <- value
312                             Matrix(x)
313                         }
314                     })
315    
316  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing",
317                                  j = "index", value = "replValue"),                                  j = "index", value = "replValue"),
318                   function(x, j, value) replDiag(x, j=j, value=value))                   function(x,i,j, ..., value) replDiag(x, j=j, value=value))
319    
320    
321  setMethod("t", signature(x = "diagonalMatrix"),  setMethod("t", signature(x = "diagonalMatrix"),
# Line 241  Line 326 
326  setMethod("isTriangular", signature(object = "diagonalMatrix"),  setMethod("isTriangular", signature(object = "diagonalMatrix"),
327            function(object) TRUE)            function(object) TRUE)
328  setMethod("isSymmetric", signature(object = "diagonalMatrix"),  setMethod("isSymmetric", signature(object = "diagonalMatrix"),
329            function(object) TRUE)            function(object, ...) TRUE)
330    
331  setMethod("chol", signature(x = "ddiMatrix"),# pivot = "ANY"  setMethod("symmpart", signature(x = "diagonalMatrix"), function(x) x)
332            function(x, pivot) {  setMethod("skewpart", signature(x = "diagonalMatrix"), setZero)
333                if(any(x@x < 0)) stop("chol() is undefined for diagonal matrix with negative entries")  
334    setMethod("chol", signature(x = "ddiMatrix"),
335              function(x, pivot, ...) {
336                  if(x@diag == "U") return(x)
337                  ## else
338                  if(any(x@x < 0))
339                      stop("chol() is undefined for diagonal matrix with negative entries")
340                x@x <- sqrt(x@x)                x@x <- sqrt(x@x)
341                x                x
342            })            })
343  ## chol(L) is L for logical diagonal:  ## chol(L) is L for logical diagonal:
344  setMethod("chol", signature(x = "ldiMatrix"), function(x, pivot) x)  setMethod("chol", signature(x = "ldiMatrix"), function(x, pivot, ...) x)
345    
346    setMethod("determinant", signature(x = "diagonalMatrix", logarithm = "logical"),
347              function(x, logarithm, ...) mkDet(.diag.x(x), logarithm))
348    
349    setMethod("norm", signature(x = "diagonalMatrix", type = "character"),
350              function(x, type, ...) {
351                  if((n <- x@Dim[1]) == 0) return(0) # as for "sparseMatrix"
352                  type <- toupper(substr(type[1], 1, 1))
353                  isU <- (x@diag == "U") # unit-diagonal
354                  if(type == "F") sqrt(if(isU) n else sum(x@x^2))
355                  else { ## norm == "I","1","O","M" :
356                      if(isU) 1 else max(abs(x@x))
357                  }
358              })
359    
360    
361    
362  ## Basic Matrix Multiplication {many more to add}  ## Basic Matrix Multiplication {many more to add}
363  ##       ---------------------  ##       ---------------------
364  ## Note that "ldi" logical are treated as numeric  ## Note that "ldi" logical are treated as numeric
365  diagdiagprod <- function(x, y) {  diagdiagprod <- function(x, y) {
366      if(any(dim(x) != dim(y))) stop("non-matching dimensions")      n <- dimCheck(x,y)[1]
367      if(x@diag != "U") {      if(x@diag != "U") {
368          if(y@diag != "U") {          if(y@diag != "U") {
369              nx <- x@x * y@x              nx <- x@x * y@x
# Line 284  Line 391 
391    
392    
393  diagmatprod <- function(x, y) {  diagmatprod <- function(x, y) {
394        ## x is diagonalMatrix
395      dx <- dim(x)      dx <- dim(x)
396      dy <- dim(y)      dy <- dim(y)
397      if(dx[2] != dy[1]) stop("non-matching dimensions")      if(dx[2] != dy[1]) stop("non-matching dimensions")
398      n <- dx[1]      n <- dx[1]
399      as(if(x@diag == "U") y else x@x * y, "Matrix")      as(if(x@diag == "U") y else x@x * y, "Matrix")
400  }  }
   
401  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),
402            diagmatprod)            diagmatprod)
403    ## sneaky .. :
404  formals(diagmatprod) <- alist(x=, y=NULL)  formals(diagmatprod) <- alist(x=, y=NULL)
405  setMethod("crossprod", signature(x = "diagonalMatrix", y = "matrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "matrix"),
406            diagmatprod)            diagmatprod)
 setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "matrix"),  
           diagmatprod)  
407    
408  diagdgeprod <- function(x, y) {  diagGeprod <- function(x, y) {
409      dx <- dim(x)      dx <- dim(x)
410      dy <- dim(y)      dy <- dim(y)
411      if(dx[2] != dy[1]) stop("non-matching dimensions")      if(dx[2] != dy[1]) stop("non-matching dimensions")
# Line 307  Line 413 
413          y@x <- x@x * y@x          y@x <- x@x * y@x
414      y      y
415  }  }
416  setMethod("%*%", signature(x = "diagonalMatrix", y = "dgeMatrix"),  setMethod("%*%", signature(x= "diagonalMatrix", y= "dgeMatrix"), diagGeprod)
417            diagdgeprod, valueClass = "dgeMatrix")  setMethod("%*%", signature(x= "diagonalMatrix", y= "lgeMatrix"), diagGeprod)
418  formals(diagdgeprod) <- alist(x=, y=NULL)  formals(diagGeprod) <- alist(x=, y=NULL)
419  setMethod("crossprod", signature(x = "diagonalMatrix", y = "dgeMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "dgeMatrix"),
420            diagdgeprod, valueClass = "dgeMatrix")            diagGeprod, valueClass = "dgeMatrix")
421    setMethod("crossprod", signature(x = "diagonalMatrix", y = "lgeMatrix"),
422              diagGeprod)
423    
424  setMethod("%*%", signature(x = "matrix", y = "diagonalMatrix"),  matdiagprod <- function(x, y) {
           function(x, y) {  
425                dx <- dim(x)                dx <- dim(x)
426                dy <- dim(y)                dy <- dim(y)
427                if(dx[2] != dy[1]) stop("non-matching dimensions")                if(dx[2] != dy[1]) stop("non-matching dimensions")
428                as(if(y@diag == "U") x else x * rep(y@x, each = dx[1]), "Matrix")      Matrix(if(y@diag == "U") x else x * rep(y@x, each = dx[1]))
429            })  }
430    setMethod("%*%", signature(x = "matrix", y = "diagonalMatrix"),
431              matdiagprod)
432    formals(matdiagprod) <- alist(x=, y=NULL)
433    setMethod("tcrossprod", signature(x = "matrix", y = "diagonalMatrix"),
434              matdiagprod)
435    
436  setMethod("%*%", signature(x = "dgeMatrix", y = "diagonalMatrix"),  gediagprod <- function(x, y) {
           function(x, y) {  
437                dx <- dim(x)                dx <- dim(x)
438                dy <- dim(y)                dy <- dim(y)
439                if(dx[2] != dy[1]) stop("non-matching dimensions")                if(dx[2] != dy[1]) stop("non-matching dimensions")
440                if(y@diag == "N")                if(y@diag == "N")
441                    x@x <- x@x * rep(y@x, each = dx[1])                    x@x <- x@x * rep(y@x, each = dx[1])
442                x                x
443            })  }
444    setMethod("%*%", signature(x= "dgeMatrix", y= "diagonalMatrix"), gediagprod)
445    setMethod("%*%", signature(x= "lgeMatrix", y= "diagonalMatrix"), gediagprod)
446    formals(gediagprod) <- alist(x=, y=NULL)
447    setMethod("tcrossprod", signature(x = "dgeMatrix", y = "diagonalMatrix"),
448              gediagprod)
449    setMethod("tcrossprod", signature(x = "lgeMatrix", y = "diagonalMatrix"),
450              gediagprod)
451    
452  ## crossprod {more of these}  ## crossprod {more of these}
453    
454  ## tcrossprod --- all are not yet there: do the dense ones here:  ## tcrossprod --- all are not yet there: do the dense ones here:
455    
456    setMethod("%*%", signature(x = "diagonalMatrix", y = "denseMatrix"),
457              function(x, y) if(x@diag == "U") y else x %*% as(y, "generalMatrix"))
458    setMethod("%*%", signature(x = "denseMatrix", y = "diagonalMatrix"),
459              function(x, y) if(y@diag == "U") x else as(x, "generalMatrix") %*% y)
460    
461    
462  ## FIXME:  ## FIXME:
463  ## setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "denseMatrix"),  ## setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "denseMatrix"),
464  ##        function(x, y = NULL) {  ##        function(x, y = NULL) {
# Line 345  Line 469 
469  ##           })  ##           })
470    
471  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
472            function(x, y = NULL) { x <- as(x, "sparseMatrix"); callGeneric() })            function(x, y = NULL) crossprod(as(x, "sparseMatrix"), y))
473    
474  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
475            function(x, y = NULL) { y <- as(y, "sparseMatrix"); callGeneric() })            function(x, y = NULL) crossprod(x, as(y, "sparseMatrix")))
476    
477  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
478            function(x, y = NULL) { x <- as(x, "sparseMatrix"); callGeneric() })            function(x, y = NULL) tcrossprod(as(x, "sparseMatrix"), y))
479    
480  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
481            function(x, y = NULL) { y <- as(y, "sparseMatrix"); callGeneric() })            function(x, y = NULL) tcrossprod(x, as(y, "sparseMatrix")))
482    
483    
484  ## 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()
485  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
486            function(x, y) as(x, "sparseMatrix") %*% y)            function(x, y) as(x, "sparseMatrix") %*% y)
487    setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
488              function(x, y) x %*% as(y, "sparseMatrix"))
489  ## NB: The previous is *not* triggering for  "ddi" o "dgC" (= distance 3)  ## NB: The previous is *not* triggering for  "ddi" o "dgC" (= distance 3)
490  ##     since there's a "ddense" o "Csparse" at dist. 2 => triggers first.  ##     since there's a "ddense" o "Csparse" at dist. 2 => triggers first.
491  ## ==> do this:  ## ==> do this:
# Line 371  Line 497 
497  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*
498  ##       do indeed work by going through sparse (and *not* ddense)!  ##       do indeed work by going through sparse (and *not* ddense)!
499    
 setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),  
           function(x, y) x %*% as(y, "sparseMatrix"))  
500    
501    
502  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),
# Line 395  Line 519 
519  setMethod("solve", signature(a = "diagonalMatrix", b = "Matrix"),  setMethod("solve", signature(a = "diagonalMatrix", b = "Matrix"),
520            solveDiag)            solveDiag)
521    
522    ## Schur()  ---> ./eigen.R
523    
524    
525    
526  ### ---------------- diagonal  o  sparse  -----------------------------  ###---------------- <Ops> (<Arith>, <Logic>, <Compare> ) ----------------------
   
527    
528  ## Use function for several signatures, in order to evade  ## Use function for several signatures, in order to evade
529  ## ambiguous dispatch for "ddi", since there's also Arith(ddense., ddense.)  ## ambiguous dispatch for "ddi", since there's also Arith(ddense., ddense.)
530  diagOdiag <- function(e1,e2) { # result should also be diagonal  diagOdiag <- function(e1,e2) {
531        ## result should also be diagonal _ if possible _
532      r <- callGeneric(.diag.x(e1), .diag.x(e2)) # error if not "compatible"      r <- callGeneric(.diag.x(e1), .diag.x(e2)) # error if not "compatible"
533        ## Check what happens with non-diagonals, i.e. (0 o 0), (FALSE o 0), ...:
534        r00 <- callGeneric(if(is.numeric(e1@x)) 0 else FALSE,
535                           if(is.numeric(e2@x)) 0 else FALSE)
536        if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* diagonal
537      if(is.numeric(r)) {      if(is.numeric(r)) {
538          if(is.numeric(e2@x)) {          if(is.numeric(e2@x)) {
539              e2@x <- r; return(.diag.2N(e2)) }              e2@x <- r; return(.diag.2N(e2)) }
# Line 418  Line 547 
547      e1@x <- r      e1@x <- r
548      .diag.2N(e1)      .diag.2N(e1)
549  }  }
550        else { ## result not diagonal, but at least symmetric:
551            isNum <- (is.numeric(r) || is.numeric(r00))
552            isLog <- (is.logical(r) || is.logical(r00))
553    
554            if(getOption("verbose"))
555                message("exploding  <diag>  o  <diag>  into dense matrix")
556            d <- e1@Dim
557            n <- d[1]
558            stopifnot(length(r) == n)
559            xx <- as.vector(matrix(rbind(r, matrix(r00,n,n)), n,n))
560            newcl <-
561                paste(if(isNum) "d" else if(isLog) {
562                    if(!any(is.na(r)) && !any(is.na(r00))) "n" else "l"
563                } else stop("not yet implemented .. please report")
564                      ,
565                      "syMatrix", sep='')
566    
567            new(newcl, Dim = e1@Dim, Dimnames = e1@Dimnames, x = xx)
568        }
569    }
570    
571    ### This would be *the* way, but we get tons of "ambiguous method dispatch"
572    ## we use this hack instead of signature  x = "diagonalMatrix" :
573    diCls <- names(getClass("diagonalMatrix")@subclasses)
574    if(FALSE) {
575  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "diagonalMatrix"),  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "diagonalMatrix"),
576            diagOdiag)            diagOdiag)
577  ## These two are just for method disambiguation:  } else { ## These are just for method disambiguation:
578  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "diagonalMatrix"),      for(c1 in diCls)
579            diagOdiag)          for(c2 in diCls)
580  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "ddiMatrix"),              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)
581            diagOdiag)  }
582    
583  ## FIXME:    diagonal  o  triangular  |-->  triangular  ## FIXME:    diagonal  o  triangular  |-->  triangular
584  ## -----     diagonal  o  symmetric   |-->  symmetric  ## -----     diagonal  o  symmetric   |-->  symmetric
# Line 435  Line 588 
588  ## For almost everything else, diag* shall be treated "as sparse" :  ## For almost everything else, diag* shall be treated "as sparse" :
589  ## These are cheap implementations via coercion  ## These are cheap implementations via coercion
590    
591  ## for disambiguation  ## For disambiguation --- define this for "sparseMatrix" , then for "ANY";
592  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "sparseMatrix"),  ## and because we can save an .M.kind() call, we use this explicit
593            function(e1,e2) callGeneric(as(e1, "sparseMatrix"), e2))  ## "hack" for all diagonalMatrix *subclasses* instead of just "diagonalMatrix" :
594  setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "diagonalMatrix"),  ##
595            function(e1,e2) callGeneric(e1, as(e2, "sparseMatrix")))  ## ddi*:
596  ## in general:  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "sparseMatrix"),
597  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "ANY"),            function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2))
598            function(e1,e2) callGeneric(as(e1,"sparseMatrix"), e2))  setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "ddiMatrix"),
599  setMethod("Ops", signature(e1 = "ANY", e2 = "diagonalMatrix"),            function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "d")))
600            function(e1,e2) callGeneric(e1, as(e2,"sparseMatrix")))  ## ldi*
601    setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "sparseMatrix"),
602              function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2))
603    setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "ldiMatrix"),
604              function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "l")))
605    
606    ## Ops:  Arith  --> numeric : "dMatrix"
607    ##       Compare --> logical
608    ##       Logic   --> logical: "lMatrix"
609    
610    ##  other = "numeric" : stay diagonal if possible
611    ## ddi*: Arith: result numeric, potentially ddiMatrix
612    setMethod("Arith", signature(e1 = "ddiMatrix", e2 = "numeric"),
613              function(e1,e2) {
614                  n <- e1@Dim[1]; nsq <- n*n
615                  f0 <- callGeneric(0, e2)
616                  if(all(is0(f0))) { # remain diagonal
617                      L1 <- (le <- length(e2)) == 1L
618                      if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)
619                      if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {
620                          e1@diag <- "N"
621                          if(L1) r <- rep.int(r, n)
622                      } else
623                          r <- callGeneric(e1@x, e2)
624                      e1@x <- if(L1) r else r[1L + n*(0:(n-1L))]
625                      return(e1)
626                  }
627                  callGeneric(diag2tT.u(e1,e2, "d"), e2)
628              })
629    
630    setMethod("Arith", signature(e1 = "numeric", e2 = "ddiMatrix"),
631              function(e1,e2) {
632                  n <- e2@Dim[1]; nsq <- n*n
633                  f0 <- callGeneric(e1, 0)
634                  if(all(is0(f0))) { # remain diagonal
635                      L1 <- (le <- length(e1)) == 1L
636                      if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)
637                      if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {
638                          e2@diag <- "N"
639                          if(L1) r <- rep.int(r, n)
640                      } else
641                          r <- callGeneric(e1, e2@x)
642                      e2@x <- if(L1) r else r[1L + n*(0:(n-1L))]
643                      return(e2)
644                  }
645                  callGeneric(e1, diag2tT.u(e2,e1, "d"))
646              })
647    
648    ## ldi* Arith --> result numeric, potentially ddiMatrix
649    setMethod("Arith", signature(e1 = "ldiMatrix", e2 = "numeric"),
650              function(e1,e2) {
651                  n <- e1@Dim[1]; nsq <- n*n
652                  f0 <- callGeneric(0, e2)
653                  if(all(is0(f0))) { # remain diagonal
654                      L1 <- (le <- length(e2)) == 1L
655                      if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)
656                      if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {
657                          e1@diag <- "N"
658                          if(L1) r <- rep.int(r, n)
659                      } else
660                          r <- callGeneric(e1@x, e2)
661                      e1 <- copyClass(e1, "ddiMatrix", c("diag", "Dim", "Dimnames"))
662                      e1@x <- if(L1) r else r[1L + n*(0:(n-1L))]
663                      return(e1)
664                  }
665                  callGeneric(diag2tT.u(e1,e2, "d"), e2)
666              })
667    
668    setMethod("Arith", signature(e1 = "numeric", e2 = "ldiMatrix"),
669              function(e1,e2) {
670                  n <- e2@Dim[1]; nsq <- n*n
671                  f0 <- callGeneric(e1, 0)
672                  if(all(is0(f0))) { # remain diagonal
673                      L1 <- (le <- length(e1)) == 1L
674                      if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)
675                      if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {
676                          e2@diag <- "N"
677                          if(L1) r <- rep.int(r, n)
678                      } else
679                          r <- callGeneric(e1, e2@x)
680                      e2 <- copyClass(e2, "ddiMatrix", c("diag", "Dim", "Dimnames"))
681                      e2@x <- if(L1) r else r[1L + n*(0:(n-1L))]
682                      return(e2)
683                  }
684                  callGeneric(e1, diag2tT.u(e2,e1, "d"))
685              })
686    
687    ## ddi*: for "Ops" without Arith --> result logical, potentially ldi
688    setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "numeric"),
689              function(e1,e2) {
690                  n <- e1@Dim[1]; nsq <- n*n
691                  f0 <- callGeneric(0, e2)
692                  if(all(is0(f0))) { # remain diagonal
693                      L1 <- (le <- length(e2)) == 1L
694                      if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)
695                      if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {
696                          e1@diag <- "N"
697                          if(L1) r <- rep.int(r, n)
698                      } else
699                          r <- callGeneric(e1@x, e2)
700                      e1 <- copyClass(e1, "ldiMatrix", c("diag", "Dim", "Dimnames"))
701                      e1@x <- if(L1) r else r[1L + n*(0:(n-1L))]
702                      return(e1)
703                  }
704                  callGeneric(diag2tT.u(e1,e2, "l"), e2)
705              })
706    
707    setMethod("Ops", signature(e1 = "numeric", e2 = "ddiMatrix"),
708              function(e1,e2) {
709                  n <- e2@Dim[1]; nsq <- n*n
710                  f0 <- callGeneric(e1, 0)
711                  if(all(is0(f0))) { # remain diagonal
712                      L1 <- (le <- length(e1)) == 1L
713                      if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)
714                      if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {
715                          e2@diag <- "N"
716                          if(L1) r <- rep.int(r, n)
717                      } else
718                          r <- callGeneric(e1, e2@x)
719                      e2 <- copyClass(e2, "ldiMatrix", c("diag", "Dim", "Dimnames"))
720                      e2@x <- if(L1) r else r[1L + n*(0:(n-1L))]
721                      return(e2)
722                  }
723                  callGeneric(e1, diag2tT.u(e2,e1, "l"))
724              })
725    
726    ## ldi*: for "Ops" without Arith --> result logical, potentially ldi
727    setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "numeric"),
728              function(e1,e2) {
729                  n <- e1@Dim[1]; nsq <- n*n
730                  f0 <- callGeneric(FALSE, e2)
731                  if(all(is0(f0))) { # remain diagonal
732                      L1 <- (le <- length(e2)) == 1L
733                      if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)
734                      if(e1@diag == "U" && any((r <- callGeneric(TRUE, e2)) != 1)) {
735                          e1@diag <- "N"
736                          if(L1) r <- rep.int(r, n)
737                      } else
738                          r <- callGeneric(e1@x, e2)
739                      e1@x <- if(L1) r else r[1L + n*(0:(n-1L))]
740                      return(e1)
741                  }
742                  callGeneric(diag2tT.u(e1,e2, "l"), e2)
743              })
744    
745    setMethod("Ops", signature(e1 = "numeric", e2 = "ldiMatrix"),
746              function(e1,e2) {
747                  n <- e2@Dim[1]; nsq <- n*n
748                  f0 <- callGeneric(e1, FALSE)
749                  if(all(is0(f0))) { # remain diagonal
750                      L1 <- (le <- length(e1)) == 1L
751                      if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)
752                      if(e2@diag == "U" && any((r <- callGeneric(e1, TRUE)) != 1)) {
753                          e2@diag <- "N"
754                          if(L1) r <- rep.int(r, n)
755                      } else
756                          r <- callGeneric(e1, e2@x)
757                      e2@x <- if(L1) r else r[1L + n*(0:(n-1L))]
758                      return(e2)
759                  }
760                  callGeneric(e1, diag2tT.u(e2,e1, "l"))
761              })
762    
763    
764    
765    ## Not {"sparseMatrix", "numeric} :  {"denseMatrix", "matrix", ... }
766    ## ddi*:
767    setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "ANY"),
768              function(e1,e2) callGeneric(diag2tT.u(e1,e2, "d"), e2))
769    setMethod("Ops", signature(e1 = "ANY", e2 = "ddiMatrix"),
770              function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "d")))
771    ## ldi*:
772    setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "ANY"),
773              function(e1,e2) callGeneric(diag2tT.u(e1,e2, "l"), e2))
774    setMethod("Ops", signature(e1 = "ANY", e2 = "ldiMatrix"),
775              function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "l")))
776    
777    ## This should *not* dispatch to <dense> methods (in ./Ops.R ), as
778    ##  FALSE & <anything> |-> FALSE : hence result should be diagonal:
779    for(cl in diCls) {
780        setMethod("&", signature(e1 = cl, e2 = "ANY"),
781                  function(e1,e2) e1 & as(e2,"Matrix"))
782        setMethod("&", signature(e1 = "ANY", e2 = cl),
783                  function(e1,e2) as(e1,"Matrix") & e2)
784        for(c2 in c("denseMatrix", "Matrix")) {
785            setMethod("&", signature(e1 = cl, e2 = c2),
786                      function(e1,e2) e1 & Diagonal(x = diag(e2)))
787            setMethod("&", signature(e1 = c2, e2 = cl),
788                      function(e1,e2) Diagonal(x = diag(e1)) & e2)
789        }
790    }
791    
792    
793    ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"
794    ### ----------  any, all: separately here
795    for(cl in diCls) {
796    setMethod("any", cl,
797              function (x, ..., na.rm) {
798                  if(any(x@Dim == 0)) FALSE
799                  else if(x@diag == "U") TRUE else any(x@x, ..., na.rm = na.rm)
800              })
801    setMethod("all",  cl, function (x, ..., na.rm) any(x@Dim == 0))
802    setMethod("prod", cl, function (x, ..., na.rm) as.numeric(any(x@Dim == 0)))
803    
804    setMethod("sum", cl,
805              function(x, ..., na.rm) {
806                  r <- sum(x@x, ..., na.rm = na.rm)# double or integer, correctly
807                  if(x@diag == "U" && !is.na(r)) r + x@Dim[1] else r
808              })
809    }
810    
811    ## The remaining ones are  max, min, range :
812    
813    setMethod("Summary", "ddiMatrix",
814              function(x, ..., na.rm) {
815                  if(any(x@Dim == 0)) callGeneric(numeric(0), ..., na.rm=na.rm)
816                  else if(x@diag == "U")
817                      callGeneric(x@x, 0, 1, ..., na.rm=na.rm)
818                  else callGeneric(x@x, 0, ..., na.rm=na.rm)
819              })
820    setMethod("Summary", "ldiMatrix",
821              function(x, ..., na.rm) {
822                  if(any(x@Dim == 0)) callGeneric(logical(0), ..., na.rm=na.rm)
823                  else if(x@diag == "U")
824                      callGeneric(x@x, FALSE, TRUE, ..., na.rm=na.rm)
825                  else callGeneric(x@x, FALSE, ..., na.rm=na.rm)
826              })
827    
828    
829    

Legend:
Removed from v.2052  
changed lines
  Added in v.2185

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