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 2096, Fri Dec 7 17:44:44 2007 UTC revision 2136, Mon Mar 17 22:20:29 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    
# Line 58  Line 60 
60      as(ret, "CsparseMatrix")      as(ret, "CsparseMatrix")
61  }  }
62    
63  diag2tT <- function(from) {  
64    .diag2tT <- function(from, uplo = "U", kind = .M.kind(from)) {
65        ## to triangular Tsparse
66      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
67      new(paste(.M.kind(from), "tTMatrix", sep=''),      new(paste(kind, "tTMatrix", sep=''),
68          diag = from@diag, Dim = from@Dim, Dimnames = from@Dimnames,          diag = from@diag, Dim = from@Dim, Dimnames = from@Dimnames,
69            uplo = uplo,
70          x = from@x, # <- ok for diag = "U" and "N" (!)          x = from@x, # <- ok for diag = "U" and "N" (!)
71          i = i, j = i)          i = i, j = i)
72  }  }
73    
74  diag2sT <- function(from) { # to symmetric Tsparse  .diag2sT <- function(from, uplo = "U", kind = .M.kind(from)) {
75      i <- if(from@diag == "U") integer(0) else seq_len(from@Dim[1]) - 1L      ## to symmetric Tsparse
76      new(paste(.M.kind(from), "sTMatrix", sep=''),      n <- from@Dim[1]
77        i <- seq_len(n) - 1L
78        new(paste(kind, "sTMatrix", sep=''),
79          Dim = from@Dim, Dimnames = from@Dimnames,          Dim = from@Dim, Dimnames = from@Dimnames,
80          x = from@x, i = i, j = i)          i = i, j = i, uplo = uplo,
81  }          x = if(from@diag == "N") from@x else ## "U"-diag
82            rep.int(switch(kind,
83  setAs("diagonalMatrix", "triangularMatrix", diag2tT)                         "d" = 1.,
84  setAs("diagonalMatrix", "sparseMatrix", diag2tT)                         "l" =,
85                           "n" = TRUE,
86                           ## otherwise
87                           stop("'", kind,"' kind not yet implemented")), n))
88    }
89    
90    ## diagonal -> triangular,  upper / lower depending on "partner":
91    diag2tT.u <- function(d, x)
92        .diag2tT(d, uplo = if(is(x,"triangularMatrix")) x@uplo else "U")
93    
94    
95    ## In order to evade method dispatch ambiguity warnings,
96    ## and because we can save a .M.kind() call, we use this explicit
97    ## "hack"  instead of signature  x = "diagonalMatrix" :
98    ##
99    ## ddi*:
100    diag2tT <- function(from) .diag2tT(from, "U", "d")
101    setAs("ddiMatrix", "triangularMatrix", diag2tT)
102    setAs("ddiMatrix", "sparseMatrix", diag2tT)
103    ## needed too (otherwise <dense> -> Tsparse is taken):
104    setAs("ddiMatrix", "TsparseMatrix", diag2tT)
105    setAs("ddiMatrix", "CsparseMatrix",
106          function(from) as(.diag2tT(from, "U", "d"), "CsparseMatrix"))
107    setAs("ddiMatrix", "symmetricMatrix",
108          function(from) .diag2sT(from, "U", "d"))
109    ##
110    ## ldi*:
111    diag2tT <- function(from) .diag2tT(from, "U", "l")
112    setAs("ldiMatrix", "triangularMatrix", diag2tT)
113    setAs("ldiMatrix", "sparseMatrix", diag2tT)
114  ## needed too (otherwise <dense> -> Tsparse is taken):  ## needed too (otherwise <dense> -> Tsparse is taken):
115  setAs("diagonalMatrix", "TsparseMatrix", diag2tT)  setAs("ldiMatrix", "TsparseMatrix", diag2tT)
116  ## is better than this:  setAs("ldiMatrix", "CsparseMatrix",
117  ## setAs("diagonalMatrix", "sparseMatrix",        function(from) as(.diag2tT(from, "U", "l"), "CsparseMatrix"))
118  ##       function(from)  setAs("ldiMatrix", "symmetricMatrix",
119  ##        as(from, if(is(from, "dMatrix")) "dgCMatrix" else "lgCMatrix"))        function(from) .diag2sT(from, "U", "l"))
120  setAs("diagonalMatrix", "CsparseMatrix",  
121        function(from) as(diag2tT(from), "CsparseMatrix"))  
122    setAs("diagonalMatrix", "nMatrix",
123          function(from) {
124              n <- from@Dim[1]
125              i <- if(from@diag == "U") integer(0) else which(isN0(from@x)) - 1L
126              new("ntTMatrix", i = i, j = i, diag = from@diag,
127                  Dim = from@Dim, Dimnames = from@Dimnames)
128          })
129    
 setAs("diagonalMatrix", "symmetricMatrix", diag2sT)  
130    
131  setAs("diagonalMatrix", "matrix",  setAs("diagonalMatrix", "matrix",
132        function(from) {        function(from) {
# Line 94  Line 136 
136                 nrow = n, ncol = n)                 nrow = n, ncol = n)
137        })        })
138    
139    setMethod("as.vector", signature(x = "diagonalMatrix", mode="missing"),
140              function(x, mode) {
141                  n <- x@Dim[1]
142                  mod <- mode(x@x)
143                  r <- vector(mod, length = n^2)
144                  if(n)
145                      r[1 + 0:(n - 1) * (n + 1)] <-
146                          if(x@diag == "U")
147                              switch(mod, "integer"= 1L,
148                                     "numeric"= 1, "logical"= TRUE)
149                          else x@x
150                  r
151              })
152    
153  setAs("diagonalMatrix", "generalMatrix", # prefer sparse:  setAs("diagonalMatrix", "generalMatrix", # prefer sparse:
154        function(from) as(from, paste(.M.kind(from), "gCMatrix", sep='')))        function(from) as(as(from, "CsparseMatrix"), "generalMatrix"))
155    
156  .diag.x <- function(m) {  .diag.x <- function(m) {
157      if(m@diag == "U")      if(m@diag == "U")
158          rep.int(if(is.numeric(m@x)) 1. else TRUE,          rep.int(if(is.numeric(m@x)) 1. else TRUE, m@Dim[1])
                 m@Dim[1])  
159      else m@x      else m@x
160  }  }
161    
# Line 109  Line 164 
164      m      m
165  }  }
166    
167    if(FALSE) {
168  ## given the above, the following  4  coercions should be all unneeded;  ## given the above, the following  4  coercions should be all unneeded;
169  ## we prefer triangular to general:  ## we prefer triangular to general:
170  setAs("ddiMatrix", "dgTMatrix",  setAs("ddiMatrix", "dgTMatrix",
# Line 139  Line 195 
195    
196  setAs("ldiMatrix", "lgCMatrix",  setAs("ldiMatrix", "lgCMatrix",
197        function(from) as(as(from, "lgTMatrix"), "lgCMatrix"))        function(from) as(as(from, "lgTMatrix"), "lgCMatrix"))
198    }
199    
200    
201  if(FALSE) # now have faster  "ddense" -> "dge"  if(FALSE) # now have faster  "ddense" -> "dge"
# Line 185  Line 242 
242        })        })
243    
244    
245  setMethod("diag", signature(x = "diagonalMatrix"),  ## In order to evade method dispatch ambiguity warnings,
246    ## we use this hack instead of signature  x = "diagonalMatrix" :
247    diCls <- names(getClass("diagonalMatrix")@subclasses)
248    for(cls in diCls) {
249        setMethod("diag", signature(x = cls),
250            function(x = 1, nrow, ncol) .diag.x(x))            function(x = 1, nrow, ncol) .diag.x(x))
251    }
252    
253    
254  subDiag <- function(x, i, j, drop) {  subDiag <- function(x, i, j, ..., drop) {
255      x <- as(x, "sparseMatrix")      x <- as(x, "sparseMatrix")
256      x <- if(missing(i))      x <- if(missing(i))
257          x[, j, drop=drop]          x[, j, drop=drop]
# Line 197  Line 259 
259          x[i, , drop=drop]          x[i, , drop=drop]
260      else      else
261          x[i,j, drop=drop]          x[i,j, drop=drop]
262      if(isDiagonal(x)) as(x, "diagonalMatrix") else x      if(isS4(x) && isDiagonal(x)) as(x, "diagonalMatrix") else x
263  }  }
264    
265  setMethod("[", signature(x = "diagonalMatrix", i = "index",  setMethod("[", signature(x = "diagonalMatrix", i = "index",
266                           j = "index", drop = "logical"), subDiag)                           j = "index", drop = "logical"), subDiag)
267  setMethod("[", signature(x = "diagonalMatrix", i = "index",  setMethod("[", signature(x = "diagonalMatrix", i = "index",
268                          j = "missing", drop = "logical"),                          j = "missing", drop = "logical"),
269            function(x, i, drop) subDiag(x, i=i, drop=drop))            function(x, i, j, ..., drop) subDiag(x, i=i, drop=drop))
270  setMethod("[", signature(x = "diagonalMatrix", i = "missing",  setMethod("[", signature(x = "diagonalMatrix", i = "missing",
271                           j = "index", drop = "logical"),                           j = "index", drop = "logical"),
272            function(x, j, drop) subDiag(x, j=j, drop=drop))            function(x, i, j, ..., drop) subDiag(x, j=j, drop=drop))
273    
274  ## When you assign to a diagonalMatrix, the result should be  ## When you assign to a diagonalMatrix, the result should be
275  ## diagonal or sparse ---  ## diagonal or sparse ---
276  ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch  ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch
277  replDiag <- function(x, i, j, value) {  ## Only(?) current bug:  x[i] <- value  is wrong when  i is *vector*
278    replDiag <- function(x, i, j, ..., value) {
279      x <- as(x, "sparseMatrix")      x <- as(x, "sparseMatrix")
     message(sprintf("diagnosing replDiag() -- nargs() = %d\n", nargs()))  
280      if(missing(i))      if(missing(i))
281          x[, j] <- value          x[, j] <- value
282      else if(missing(j))      else if(missing(j)) { ##  x[i , ] <- v  *OR*   x[i] <- v
283            na <- nargs()
284    ##         message("diagnosing replDiag() -- nargs()= ", na)
285            if(na == 4)
286          x[i, ] <- value          x[i, ] <- value
287      else          else if(na == 3)
288                x[i] <- value
289            else stop("Internal bug: nargs()=",na,"; please report")
290        } else
291          x[i,j] <- value          x[i,j] <- value
292      if(isDiagonal(x)) as(x, "diagonalMatrix") else x      if(isDiagonal(x)) as(x, "diagonalMatrix") else x
293  }  }
294    
295  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index",
296                                  j = "index", value = "replValue"), replDiag)                                  j = "index", value = "replValue"), replDiag)
297    
298  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index",
299                                  j = "missing", value = "replValue"),                                  j = "missing", value = "replValue"),
300                   function(x, i, value) replDiag(x, i=i, value=value))                   function(x,i,j, ..., value) {
301                         ## message("before replDiag() -- nargs()= ", nargs())
302                         if(nargs() == 3)
303                             replDiag(x, i=i, value=value)
304                         else ## nargs() == 4 :
305                             replDiag(x, i=i, , value=value)
306                     })
307    
308  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "matrix", # 2-col.matrix  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "matrix", # 2-col.matrix
309                                  j = "missing", value = "replValue"),                                  j = "missing", value = "replValue"),
310                   function(x, i, value) {                   function(x,i,j, ..., value) {
311                       if(ncol(i) == 2) {                       if(ncol(i) == 2) {
312                           if(all((ii <- i[,1]) == i[,2])) { # replace in diagonal only                           if(all((ii <- i[,1]) == i[,2])) { # replace in diagonal only
313                               x@x[ii] <- value                               x@x[ii] <- value
# Line 251  Line 327 
327    
328  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing",
329                                  j = "index", value = "replValue"),                                  j = "index", value = "replValue"),
330                   function(x, j, value) replDiag(x, j=j, value=value))                   function(x,i,j, ..., value) replDiag(x, j=j, value=value))
331    
332    
333  setMethod("t", signature(x = "diagonalMatrix"),  setMethod("t", signature(x = "diagonalMatrix"),
# Line 262  Line 338 
338  setMethod("isTriangular", signature(object = "diagonalMatrix"),  setMethod("isTriangular", signature(object = "diagonalMatrix"),
339            function(object) TRUE)            function(object) TRUE)
340  setMethod("isSymmetric", signature(object = "diagonalMatrix"),  setMethod("isSymmetric", signature(object = "diagonalMatrix"),
341            function(object) TRUE)            function(object, ...) TRUE)
342    
343    setMethod("symmpart", signature(x = "diagonalMatrix"), function(x) x)
344    setMethod("skewpart", signature(x = "diagonalMatrix"), setZero)
345    
346  setMethod("chol", signature(x = "ddiMatrix"),# pivot = "ANY"  setMethod("chol", signature(x = "ddiMatrix"),
347            function(x, pivot) {            function(x, pivot, ...) {
348                if(any(x@x < 0)) stop("chol() is undefined for diagonal matrix with negative entries")                if(any(x@x < 0))
349                      stop("chol() is undefined for diagonal matrix with negative entries")
350                x@x <- sqrt(x@x)                x@x <- sqrt(x@x)
351                x                x
352            })            })
353  ## chol(L) is L for logical diagonal:  ## chol(L) is L for logical diagonal:
354  setMethod("chol", signature(x = "ldiMatrix"), function(x, pivot) x)  setMethod("chol", signature(x = "ldiMatrix"), function(x, pivot, ...) x)
355    
356  ## Basic Matrix Multiplication {many more to add}  ## Basic Matrix Multiplication {many more to add}
357  ##       ---------------------  ##       ---------------------
# Line 366  Line 446 
446  ##           })  ##           })
447    
448  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
449            function(x, y = NULL) { x <- as(x, "sparseMatrix"); callGeneric() })            function(x, y = NULL) crossprod(as(x, "sparseMatrix"), y))
450    
451  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
452            function(x, y = NULL) { y <- as(y, "sparseMatrix"); callGeneric() })            function(x, y = NULL) crossprod(x, as(y, "sparseMatrix")))
453    
454  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
455            function(x, y = NULL) { x <- as(x, "sparseMatrix"); callGeneric() })            function(x, y = NULL) tcrossprod(as(x, "sparseMatrix"), y))
456    
457  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
458            function(x, y = NULL) { y <- as(y, "sparseMatrix"); callGeneric() })            function(x, y = NULL) tcrossprod(x, as(y, "sparseMatrix")))
459    
460    
461  ## 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()
# Line 416  Line 496 
496  setMethod("solve", signature(a = "diagonalMatrix", b = "Matrix"),  setMethod("solve", signature(a = "diagonalMatrix", b = "Matrix"),
497            solveDiag)            solveDiag)
498    
499    ## Schur()  ---> ./eigen.R
500    
501    
502    
# Line 456  Line 537 
537  ## For almost everything else, diag* shall be treated "as sparse" :  ## For almost everything else, diag* shall be treated "as sparse" :
538  ## These are cheap implementations via coercion  ## These are cheap implementations via coercion
539    
540  ## for disambiguation  ## for disambiguation --- define this for "sparseMatrix" , then for "ANY" :
541  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "sparseMatrix"),  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "sparseMatrix"),
542            function(e1,e2) callGeneric(as(e1, "sparseMatrix"), e2))            function(e1,e2) callGeneric(diag2tT.u(e1,e2), e2))
543  setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "diagonalMatrix"),  setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "diagonalMatrix"),
544            function(e1,e2) callGeneric(e1, as(e2, "sparseMatrix")))            function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1)))
545  ## in general:  ## in general:
546  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "ANY"),  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "ANY"),
547            function(e1,e2) callGeneric(as(e1,"sparseMatrix"), e2))            function(e1,e2) callGeneric(diag2tT.u(e1,e2), e2))
548  setMethod("Ops", signature(e1 = "ANY", e2 = "diagonalMatrix"),  setMethod("Ops", signature(e1 = "ANY", e2 = "diagonalMatrix"),
549            function(e1,e2) callGeneric(e1, as(e2,"sparseMatrix")))            function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1)))
550    
551    
552    

Legend:
Removed from v.2096  
changed lines
  Added in v.2136

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