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 2203, Sat Jun 14 20:09:17 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",
297                                    i = "matrix", # 2-col.matrix
298                                    j = "missing", value = "replValue"),
299                     function(x,i,j, ..., value) {
300                         if(ncol(i) == 2) {
301                             if(all((ii <- i[,1]) == i[,2])) { # replace in diagonal only
302                                 if(x@diag == "U") {
303                                     one <- as1(x@x)
304                                     if(any(value != one | is.na(value))) {
305                                         x@diag <- "N"
306                                         x@x <- rep.int(one, x@Dim[1])
307                                     }
308                                 }
309                                 x@x[ii] <- value
310                                 x
311                             } else { ## no longer diagonal, but remain sparse:
312                                 x <- as(x, "sparseMatrix")
313                                 x[i] <- value
314                                 x
315                             }
316                         }
317                         else { # behave as "base R": use as if vector
318                             x <- as(x, "matrix")
319                             x[i] <- value
320                             Matrix(x)
321                         }
322                     })
323    
324  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing",
325                                  j = "index", value = "replValue"),                                  j = "index", value = "replValue"),
326                   function(x, j, value) replDiag(x, j=j, value=value))                   function(x,i,j, ..., value) replDiag(x, j=j, value=value))
327    
328    
329  setMethod("t", signature(x = "diagonalMatrix"),  setMethod("t", signature(x = "diagonalMatrix"),
# Line 241  Line 334 
334  setMethod("isTriangular", signature(object = "diagonalMatrix"),  setMethod("isTriangular", signature(object = "diagonalMatrix"),
335            function(object) TRUE)            function(object) TRUE)
336  setMethod("isSymmetric", signature(object = "diagonalMatrix"),  setMethod("isSymmetric", signature(object = "diagonalMatrix"),
337            function(object) TRUE)            function(object, ...) TRUE)
338    
339  setMethod("chol", signature(x = "ddiMatrix"),# pivot = "ANY"  setMethod("symmpart", signature(x = "diagonalMatrix"), function(x) x)
340            function(x, pivot) {  setMethod("skewpart", signature(x = "diagonalMatrix"), setZero)
341                if(any(x@x < 0)) stop("chol() is undefined for diagonal matrix with negative entries")  
342    setMethod("chol", signature(x = "ddiMatrix"),
343              function(x, pivot, ...) {
344                  if(x@diag == "U") return(x)
345                  ## else
346                  if(any(x@x < 0))
347                      stop("chol() is undefined for diagonal matrix with negative entries")
348                x@x <- sqrt(x@x)                x@x <- sqrt(x@x)
349                x                x
350            })            })
351  ## chol(L) is L for logical diagonal:  ## chol(L) is L for logical diagonal:
352  setMethod("chol", signature(x = "ldiMatrix"), function(x, pivot) x)  setMethod("chol", signature(x = "ldiMatrix"), function(x, pivot, ...) x)
353    
354    setMethod("determinant", signature(x = "diagonalMatrix", logarithm = "logical"),
355              function(x, logarithm, ...) mkDet(.diag.x(x), logarithm))
356    
357    setMethod("norm", signature(x = "diagonalMatrix", type = "character"),
358              function(x, type, ...) {
359                  if((n <- x@Dim[1]) == 0) return(0) # as for "sparseMatrix"
360                  type <- toupper(substr(type[1], 1, 1))
361                  isU <- (x@diag == "U") # unit-diagonal
362                  if(type == "F") sqrt(if(isU) n else sum(x@x^2))
363                  else { ## norm == "I","1","O","M" :
364                      if(isU) 1 else max(abs(x@x))
365                  }
366              })
367    
368    
369    
370  ## Basic Matrix Multiplication {many more to add}  ## Basic Matrix Multiplication {many more to add}
371  ##       ---------------------  ##       ---------------------
372  ## Note that "ldi" logical are treated as numeric  ## Note that "ldi" logical are treated as numeric
373  diagdiagprod <- function(x, y) {  diagdiagprod <- function(x, y) {
374      if(any(dim(x) != dim(y))) stop("non-matching dimensions")      n <- dimCheck(x,y)[1]
375      if(x@diag != "U") {      if(x@diag != "U") {
376          if(y@diag != "U") {          if(y@diag != "U") {
377              nx <- x@x * y@x              nx <- x@x * y@x
# Line 284  Line 399 
399    
400    
401  diagmatprod <- function(x, y) {  diagmatprod <- function(x, y) {
402        ## x is diagonalMatrix
403      dx <- dim(x)      dx <- dim(x)
404      dy <- dim(y)      dy <- dim(y)
405      if(dx[2] != dy[1]) stop("non-matching dimensions")      if(dx[2] != dy[1]) stop("non-matching dimensions")
406      n <- dx[1]      n <- dx[1]
407      as(if(x@diag == "U") y else x@x * y, "Matrix")      as(if(x@diag == "U") y else x@x * y, "Matrix")
408  }  }
   
409  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),
410            diagmatprod)            diagmatprod)
411    ## sneaky .. :
412  formals(diagmatprod) <- alist(x=, y=NULL)  formals(diagmatprod) <- alist(x=, y=NULL)
413  setMethod("crossprod", signature(x = "diagonalMatrix", y = "matrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "matrix"),
414            diagmatprod)            diagmatprod)
 setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "matrix"),  
           diagmatprod)  
415    
416  diagdgeprod <- function(x, y) {  diagGeprod <- function(x, y) {
417      dx <- dim(x)      dx <- dim(x)
418      dy <- dim(y)      dy <- dim(y)
419      if(dx[2] != dy[1]) stop("non-matching dimensions")      if(dx[2] != dy[1]) stop("non-matching dimensions")
# Line 307  Line 421 
421          y@x <- x@x * y@x          y@x <- x@x * y@x
422      y      y
423  }  }
424  setMethod("%*%", signature(x = "diagonalMatrix", y = "dgeMatrix"),  setMethod("%*%", signature(x= "diagonalMatrix", y= "dgeMatrix"), diagGeprod)
425            diagdgeprod, valueClass = "dgeMatrix")  setMethod("%*%", signature(x= "diagonalMatrix", y= "lgeMatrix"), diagGeprod)
426  formals(diagdgeprod) <- alist(x=, y=NULL)  formals(diagGeprod) <- alist(x=, y=NULL)
427  setMethod("crossprod", signature(x = "diagonalMatrix", y = "dgeMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "dgeMatrix"),
428            diagdgeprod, valueClass = "dgeMatrix")            diagGeprod, valueClass = "dgeMatrix")
429    setMethod("crossprod", signature(x = "diagonalMatrix", y = "lgeMatrix"),
430              diagGeprod)
431    
432  setMethod("%*%", signature(x = "matrix", y = "diagonalMatrix"),  matdiagprod <- function(x, y) {
           function(x, y) {  
433                dx <- dim(x)                dx <- dim(x)
434                dy <- dim(y)                dy <- dim(y)
435                if(dx[2] != dy[1]) stop("non-matching dimensions")                if(dx[2] != dy[1]) stop("non-matching dimensions")
436                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]))
437            })  }
438    setMethod("%*%", signature(x = "matrix", y = "diagonalMatrix"),
439              matdiagprod)
440    formals(matdiagprod) <- alist(x=, y=NULL)
441    setMethod("tcrossprod", signature(x = "matrix", y = "diagonalMatrix"),
442              matdiagprod)
443    
444  setMethod("%*%", signature(x = "dgeMatrix", y = "diagonalMatrix"),  gediagprod <- function(x, y) {
           function(x, y) {  
445                dx <- dim(x)                dx <- dim(x)
446                dy <- dim(y)                dy <- dim(y)
447                if(dx[2] != dy[1]) stop("non-matching dimensions")                if(dx[2] != dy[1]) stop("non-matching dimensions")
448                if(y@diag == "N")                if(y@diag == "N")
449                    x@x <- x@x * rep(y@x, each = dx[1])                    x@x <- x@x * rep(y@x, each = dx[1])
450                x                x
451            })  }
452    setMethod("%*%", signature(x= "dgeMatrix", y= "diagonalMatrix"), gediagprod)
453    setMethod("%*%", signature(x= "lgeMatrix", y= "diagonalMatrix"), gediagprod)
454    formals(gediagprod) <- alist(x=, y=NULL)
455    setMethod("tcrossprod", signature(x = "dgeMatrix", y = "diagonalMatrix"),
456              gediagprod)
457    setMethod("tcrossprod", signature(x = "lgeMatrix", y = "diagonalMatrix"),
458              gediagprod)
459    
460  ## crossprod {more of these}  ## crossprod {more of these}
461    
462  ## tcrossprod --- all are not yet there: do the dense ones here:  ## tcrossprod --- all are not yet there: do the dense ones here:
463    
464    setMethod("%*%", signature(x = "diagonalMatrix", y = "denseMatrix"),
465              function(x, y) if(x@diag == "U") y else x %*% as(y, "generalMatrix"))
466    setMethod("%*%", signature(x = "denseMatrix", y = "diagonalMatrix"),
467              function(x, y) if(y@diag == "U") x else as(x, "generalMatrix") %*% y)
468    
469    
470  ## FIXME:  ## FIXME:
471  ## setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "denseMatrix"),  ## setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "denseMatrix"),
472  ##        function(x, y = NULL) {  ##        function(x, y = NULL) {
# Line 345  Line 477 
477  ##           })  ##           })
478    
479  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
480            function(x, y = NULL) { x <- as(x, "sparseMatrix"); callGeneric() })            function(x, y = NULL) crossprod(as(x, "sparseMatrix"), y))
481    
482  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
483            function(x, y = NULL) { y <- as(y, "sparseMatrix"); callGeneric() })            function(x, y = NULL) crossprod(x, as(y, "sparseMatrix")))
484    
485  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
486            function(x, y = NULL) { x <- as(x, "sparseMatrix"); callGeneric() })            function(x, y = NULL) tcrossprod(as(x, "sparseMatrix"), y))
487    
488  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
489            function(x, y = NULL) { y <- as(y, "sparseMatrix"); callGeneric() })            function(x, y = NULL) tcrossprod(x, as(y, "sparseMatrix")))
490    
491    
492  ## 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()
493  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
494            function(x, y) as(x, "sparseMatrix") %*% y)            function(x, y) as(x, "sparseMatrix") %*% y)
495    setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
496              function(x, y) x %*% as(y, "sparseMatrix"))
497  ## NB: The previous is *not* triggering for  "ddi" o "dgC" (= distance 3)  ## NB: The previous is *not* triggering for  "ddi" o "dgC" (= distance 3)
498  ##     since there's a "ddense" o "Csparse" at dist. 2 => triggers first.  ##     since there's a "ddense" o "Csparse" at dist. 2 => triggers first.
499  ## ==> do this:  ## ==> do this:
# Line 371  Line 505 
505  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*
506  ##       do indeed work by going through sparse (and *not* ddense)!  ##       do indeed work by going through sparse (and *not* ddense)!
507    
 setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),  
           function(x, y) x %*% as(y, "sparseMatrix"))  
508    
509    
510  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),
# Line 395  Line 527 
527  setMethod("solve", signature(a = "diagonalMatrix", b = "Matrix"),  setMethod("solve", signature(a = "diagonalMatrix", b = "Matrix"),
528            solveDiag)            solveDiag)
529    
530    ## Schur()  ---> ./eigen.R
531    
532    
533    
534  ### ---------------- diagonal  o  sparse  -----------------------------  ###---------------- <Ops> (<Arith>, <Logic>, <Compare> ) ----------------------
   
535    
536  ## Use function for several signatures, in order to evade  ## Use function for several signatures, in order to evade
537  ## ambiguous dispatch for "ddi", since there's also Arith(ddense., ddense.)  ## ambiguous dispatch for "ddi", since there's also Arith(ddense., ddense.)
538  diagOdiag <- function(e1,e2) { # result should also be diagonal  diagOdiag <- function(e1,e2) {
539        ## result should also be diagonal _ if possible _
540      r <- callGeneric(.diag.x(e1), .diag.x(e2)) # error if not "compatible"      r <- callGeneric(.diag.x(e1), .diag.x(e2)) # error if not "compatible"
541        ## Check what happens with non-diagonals, i.e. (0 o 0), (FALSE o 0), ...:
542        r00 <- callGeneric(if(is.numeric(e1@x)) 0 else FALSE,
543                           if(is.numeric(e2@x)) 0 else FALSE)
544        if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* diagonal
545      if(is.numeric(r)) {      if(is.numeric(r)) {
546          if(is.numeric(e2@x)) {          if(is.numeric(e2@x)) {
547              e2@x <- r; return(.diag.2N(e2)) }              e2@x <- r; return(.diag.2N(e2)) }
# Line 418  Line 555 
555      e1@x <- r      e1@x <- r
556      .diag.2N(e1)      .diag.2N(e1)
557  }  }
558        else { ## result not diagonal, but at least symmetric:
559            isNum <- (is.numeric(r) || is.numeric(r00))
560            isLog <- (is.logical(r) || is.logical(r00))
561    
562            if(getOption("verbose"))
563                message("exploding  <diag>  o  <diag>  into dense matrix")
564            d <- e1@Dim
565            n <- d[1]
566            stopifnot(length(r) == n)
567            xx <- as.vector(matrix(rbind(r, matrix(r00,n,n)), n,n))
568            newcl <-
569                paste(if(isNum) "d" else if(isLog) {
570                    if(!any(is.na(r)) && !any(is.na(r00))) "n" else "l"
571                } else stop("not yet implemented .. please report")
572                      ,
573                      "syMatrix", sep='')
574    
575            new(newcl, Dim = e1@Dim, Dimnames = e1@Dimnames, x = xx)
576        }
577    }
578    
579    ### This would be *the* way, but we get tons of "ambiguous method dispatch"
580    ## we use this hack instead of signature  x = "diagonalMatrix" :
581    diCls <- names(getClass("diagonalMatrix")@subclasses)
582    if(FALSE) {
583  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "diagonalMatrix"),  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "diagonalMatrix"),
584            diagOdiag)            diagOdiag)
585  ## These two are just for method disambiguation:  } else { ## These are just for method disambiguation:
586  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "diagonalMatrix"),      for(c1 in diCls)
587            diagOdiag)          for(c2 in diCls)
588  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "ddiMatrix"),              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)
589            diagOdiag)  }
590    
591  ## FIXME:    diagonal  o  triangular  |-->  triangular  ## FIXME:    diagonal  o  triangular  |-->  triangular
592  ## -----     diagonal  o  symmetric   |-->  symmetric  ## -----     diagonal  o  symmetric   |-->  symmetric
# Line 435  Line 596 
596  ## For almost everything else, diag* shall be treated "as sparse" :  ## For almost everything else, diag* shall be treated "as sparse" :
597  ## These are cheap implementations via coercion  ## These are cheap implementations via coercion
598    
599  ## for disambiguation  ## For disambiguation --- define this for "sparseMatrix" , then for "ANY";
600  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "sparseMatrix"),  ## and because we can save an .M.kind() call, we use this explicit
601            function(e1,e2) callGeneric(as(e1, "sparseMatrix"), e2))  ## "hack" for all diagonalMatrix *subclasses* instead of just "diagonalMatrix" :
602  setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "diagonalMatrix"),  ##
603            function(e1,e2) callGeneric(e1, as(e2, "sparseMatrix")))  ## ddi*:
604  ## in general:  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "sparseMatrix"),
605  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "ANY"),            function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2))
606            function(e1,e2) callGeneric(as(e1,"sparseMatrix"), e2))  setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "ddiMatrix"),
607  setMethod("Ops", signature(e1 = "ANY", e2 = "diagonalMatrix"),            function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "d")))
608            function(e1,e2) callGeneric(e1, as(e2,"sparseMatrix")))  ## ldi*
609    setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "sparseMatrix"),
610              function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2))
611    setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "ldiMatrix"),
612              function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "l")))
613    
614    ## Ops:  Arith  --> numeric : "dMatrix"
615    ##       Compare --> logical
616    ##       Logic   --> logical: "lMatrix"
617    
618    ##  other = "numeric" : stay diagonal if possible
619    ## ddi*: Arith: result numeric, potentially ddiMatrix
620    setMethod("Arith", signature(e1 = "ddiMatrix", e2 = "numeric"),
621              function(e1,e2) {
622                  n <- e1@Dim[1]; nsq <- n^2
623                  f0 <- callGeneric(0, e2)
624                  if(all(is0(f0))) { # remain diagonal
625                      L1 <- (le <- length(e2)) == 1L
626                      if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)
627                      if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {
628                          e1@diag <- "N"
629                          if(L1) r <- rep.int(r, n)
630                      } else
631                          r <- callGeneric(e1@x, e2)
632                      e1@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
633                      return(e1)
634                  }
635                  callGeneric(diag2tT.u(e1,e2, "d"), e2)
636              })
637    
638    setMethod("Arith", signature(e1 = "numeric", e2 = "ddiMatrix"),
639              function(e1,e2) {
640                  n <- e2@Dim[1]; nsq <- n^2
641                  f0 <- callGeneric(e1, 0)
642                  if(all(is0(f0))) { # remain diagonal
643                      L1 <- (le <- length(e1)) == 1L
644                      if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)
645                      if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {
646                          e2@diag <- "N"
647                          if(L1) r <- rep.int(r, n)
648                      } else
649                          r <- callGeneric(e1, e2@x)
650                      e2@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
651                      return(e2)
652                  }
653                  callGeneric(e1, diag2tT.u(e2,e1, "d"))
654              })
655    
656    ## ldi* Arith --> result numeric, potentially ddiMatrix
657    setMethod("Arith", signature(e1 = "ldiMatrix", e2 = "numeric"),
658              function(e1,e2) {
659                  n <- e1@Dim[1]; nsq <- n^2
660                  f0 <- callGeneric(0, e2)
661                  if(all(is0(f0))) { # remain diagonal
662                      L1 <- (le <- length(e2)) == 1L
663                      if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)
664                      if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {
665                          e1@diag <- "N"
666                          if(L1) r <- rep.int(r, n)
667                      } else
668                          r <- callGeneric(e1@x, e2)
669                      e1 <- copyClass(e1, "ddiMatrix", c("diag", "Dim", "Dimnames"))
670                      e1@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
671                      return(e1)
672                  }
673                  callGeneric(diag2tT.u(e1,e2, "d"), e2)
674              })
675    
676    setMethod("Arith", signature(e1 = "numeric", e2 = "ldiMatrix"),
677              function(e1,e2) {
678                  n <- e2@Dim[1]; nsq <- n^2
679                  f0 <- callGeneric(e1, 0)
680                  if(all(is0(f0))) { # remain diagonal
681                      L1 <- (le <- length(e1)) == 1L
682                      if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)
683                      if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {
684                          e2@diag <- "N"
685                          if(L1) r <- rep.int(r, n)
686                      } else
687                          r <- callGeneric(e1, e2@x)
688                      e2 <- copyClass(e2, "ddiMatrix", c("diag", "Dim", "Dimnames"))
689                      e2@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
690                      return(e2)
691                  }
692                  callGeneric(e1, diag2tT.u(e2,e1, "d"))
693              })
694    
695    ## ddi*: for "Ops" without Arith --> result logical, potentially ldi
696    setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "numeric"),
697              function(e1,e2) {
698                  n <- e1@Dim[1]; nsq <- n^2
699                  f0 <- callGeneric(0, e2)
700                  if(all(is0(f0))) { # remain diagonal
701                      L1 <- (le <- length(e2)) == 1L
702                      if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)
703                      if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {
704                          e1@diag <- "N"
705                          if(L1) r <- rep.int(r, n)
706                      } else
707                          r <- callGeneric(e1@x, e2)
708                      e1 <- copyClass(e1, "ldiMatrix", c("diag", "Dim", "Dimnames"))
709                      e1@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
710                      return(e1)
711                  }
712                  callGeneric(diag2tT.u(e1,e2, "l"), e2)
713              })
714    
715    setMethod("Ops", signature(e1 = "numeric", e2 = "ddiMatrix"),
716              function(e1,e2) {
717                  n <- e2@Dim[1]; nsq <- n^2
718                  f0 <- callGeneric(e1, 0)
719                  if(all(is0(f0))) { # remain diagonal
720                      L1 <- (le <- length(e1)) == 1L
721                      if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)
722                      if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {
723                          e2@diag <- "N"
724                          if(L1) r <- rep.int(r, n)
725                      } else
726                          r <- callGeneric(e1, e2@x)
727                      e2 <- copyClass(e2, "ldiMatrix", c("diag", "Dim", "Dimnames"))
728                      e2@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
729                      return(e2)
730                  }
731                  callGeneric(e1, diag2tT.u(e2,e1, "l"))
732              })
733    
734    ## ldi*: for "Ops" without Arith --> result logical, potentially ldi
735    setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "numeric"),
736              function(e1,e2) {
737                  n <- e1@Dim[1]; nsq <- n^2
738                  f0 <- callGeneric(FALSE, e2)
739                  if(all(is0(f0))) { # remain diagonal
740                      L1 <- (le <- length(e2)) == 1L
741                      if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)
742                      if(e1@diag == "U" && any((r <- callGeneric(TRUE, e2)) != 1)) {
743                          e1@diag <- "N"
744                          if(L1) r <- rep.int(r, n)
745                      } else
746                          r <- callGeneric(e1@x, e2)
747                      e1@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
748                      return(e1)
749                  }
750                  callGeneric(diag2tT.u(e1,e2, "l"), e2)
751              })
752    
753    setMethod("Ops", signature(e1 = "numeric", e2 = "ldiMatrix"),
754              function(e1,e2) {
755                  n <- e2@Dim[1]; nsq <- n^2
756                  f0 <- callGeneric(e1, FALSE)
757                  if(all(is0(f0))) { # remain diagonal
758                      L1 <- (le <- length(e1)) == 1L
759                      if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)
760                      if(e2@diag == "U" && any((r <- callGeneric(e1, TRUE)) != 1)) {
761                          e2@diag <- "N"
762                          if(L1) r <- rep.int(r, n)
763                      } else
764                          r <- callGeneric(e1, e2@x)
765                      e2@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
766                      return(e2)
767                  }
768                  callGeneric(e1, diag2tT.u(e2,e1, "l"))
769              })
770    
771    
772    
773    ## Not {"sparseMatrix", "numeric} :  {"denseMatrix", "matrix", ... }
774    for(other in c("ANY", "Matrix", "dMatrix")) {
775        ## ddi*:
776        setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other),
777                  function(e1,e2) callGeneric(diag2tT.u(e1,e2, "d"), e2))
778        setMethod("Ops", signature(e1 = other, e2 = "ddiMatrix"),
779                  function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "d")))
780        ## ldi*:
781        setMethod("Ops", signature(e1 = "ldiMatrix", e2 = other),
782                  function(e1,e2) callGeneric(diag2tT.u(e1,e2, "l"), e2))
783        setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),
784                  function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "l")))
785    }
786    ## This should *not* dispatch to <dense> methods (in ./Ops.R ), as
787    ##  FALSE & <anything> |-> FALSE : hence result should be diagonal:
788    for(cl in diCls) {
789        setMethod("&", signature(e1 = cl, e2 = "ANY"),
790                  function(e1,e2) e1 & as(e2,"Matrix"))
791        setMethod("&", signature(e1 = "ANY", e2 = cl),
792                  function(e1,e2) as(e1,"Matrix") & e2)
793        for(c2 in c("denseMatrix", "Matrix")) {
794            setMethod("&", signature(e1 = cl, e2 = c2),
795                      function(e1,e2) e1 & Diagonal(x = diag(e2)))
796            setMethod("&", signature(e1 = c2, e2 = cl),
797                      function(e1,e2) Diagonal(x = diag(e1)) & e2)
798        }
799    }
800    
801    
802    ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"
803    ### ----------  any, all: separately here
804    for(cl in diCls) {
805    setMethod("any", cl,
806              function (x, ..., na.rm) {
807                  if(any(x@Dim == 0)) FALSE
808                  else if(x@diag == "U") TRUE else any(x@x, ..., na.rm = na.rm)
809              })
810    setMethod("all",  cl, function (x, ..., na.rm) any(x@Dim == 0))
811    setMethod("prod", cl, function (x, ..., na.rm) as.numeric(any(x@Dim == 0)))
812    
813    setMethod("sum", cl,
814              function(x, ..., na.rm) {
815                  r <- sum(x@x, ..., na.rm = na.rm)# double or integer, correctly
816                  if(x@diag == "U" && !is.na(r)) r + x@Dim[1] else r
817              })
818    }
819    
820    ## The remaining ones are  max, min, range :
821    
822    setMethod("Summary", "ddiMatrix",
823              function(x, ..., na.rm) {
824                  if(any(x@Dim == 0)) callGeneric(numeric(0), ..., na.rm=na.rm)
825                  else if(x@diag == "U")
826                      callGeneric(x@x, 0, 1, ..., na.rm=na.rm)
827                  else callGeneric(x@x, 0, ..., na.rm=na.rm)
828              })
829    setMethod("Summary", "ldiMatrix",
830              function(x, ..., na.rm) {
831                  if(any(x@Dim == 0)) callGeneric(logical(0), ..., na.rm=na.rm)
832                  else if(x@diag == "U")
833                      callGeneric(x@x, FALSE, TRUE, ..., na.rm=na.rm)
834                  else callGeneric(x@x, FALSE, ..., na.rm=na.rm)
835              })
836    
837    
838    
# Line 463  Line 851 
851            function(object) {            function(object) {
852                d <- dim(object)                d <- dim(object)
853                cl <- class(object)                cl <- class(object)
854                cat(sprintf('%d x %d diagonal matrix of class "%s"\n',                cat(sprintf('%d x %d diagonal matrix of class "%s"',
855                            d[1], d[2], cl))                            d[1], d[2], cl))
856                  if(d[1] < 50) {
857                      cat("\n")
858                prDiag(object)                prDiag(object)
859                  } else {
860                      cat(", with diagonal entries\n")
861                      show(diag(object))
862                      invisible(object)
863                  }
864            })            })

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

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