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 2128, Thu Mar 13 23:08:34 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") { ## to triangular Tsparse
65      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
66      new(paste(.M.kind(from), "tTMatrix", sep=''),      new(paste(.M.kind(from), "tTMatrix", sep=''),
67          diag = from@diag, Dim = from@Dim, Dimnames = from@Dimnames,          diag = from@diag, Dim = from@Dim, Dimnames = from@Dimnames,
68            uplo = uplo,
69          x = from@x, # <- ok for diag = "U" and "N" (!)          x = from@x, # <- ok for diag = "U" and "N" (!)
70          i = i, j = i)          i = i, j = i)
71  }  }
72    
73  diag2sT <- function(from) { # to symmetric Tsparse  diag2tT <- function(from) .diag2tT(from, "U")
74      i <- if(from@diag == "U") integer(0) else seq_len(from@Dim[1]) - 1L  
75      new(paste(.M.kind(from), "sTMatrix", sep=''),  .diag2sT <- function(from, uplo = "U") { # to symmetric Tsparse
76        n <- from@Dim[1]
77        i <- seq_len(n) - 1L
78        kind <- .M.kind(from)
79        new(paste(kind, "sTMatrix", sep=''),
80          Dim = from@Dim, Dimnames = from@Dimnames,          Dim = from@Dim, Dimnames = from@Dimnames,
81          x = from@x, i = i, j = i)          i = i, j = i, uplo = uplo,
82            x = if(from@diag == "N") from@x else ## "U"-diag
83            rep.int(switch(kind,
84                           "d" = 1.,
85                           "l" =,
86                           "n" = TRUE,
87                           ## otherwise
88                           stop("'", kind,"' kind not yet implemented")), n))
89  }  }
90    
91    diag2sT <- function(from) .diag2sT(from, "U")
92    
93    ## diagonal -> triangular,  upper / lower depending on "partner":
94    diag2tT.u <- function(d, x)
95        .diag2tT(d, uplo = if(is(x,"triangularMatrix")) x@uplo else "U")
96    
97  setAs("diagonalMatrix", "triangularMatrix", diag2tT)  setAs("diagonalMatrix", "triangularMatrix", diag2tT)
98  setAs("diagonalMatrix", "sparseMatrix", diag2tT)  setAs("diagonalMatrix", "sparseMatrix", diag2tT)
99  ## needed too (otherwise <dense> -> Tsparse is taken):  ## needed too (otherwise <dense> -> Tsparse is taken):
# Line 86  Line 107 
107    
108  setAs("diagonalMatrix", "symmetricMatrix", diag2sT)  setAs("diagonalMatrix", "symmetricMatrix", diag2sT)
109    
110    setAs("diagonalMatrix", "nMatrix",
111          function(from) {
112              n <- from@Dim[1]
113              i <- if(from@diag == "U") integer(0) else which(isN0(from@x)) - 1L
114              new("ntTMatrix", i = i, j = i, diag = from@diag,
115                  Dim = from@Dim, Dimnames = from@Dimnames)
116          })
117    
118    
119  setAs("diagonalMatrix", "matrix",  setAs("diagonalMatrix", "matrix",
120        function(from) {        function(from) {
121            n <- from@Dim[1]            n <- from@Dim[1]
# Line 94  Line 124 
124                 nrow = n, ncol = n)                 nrow = n, ncol = n)
125        })        })
126    
127    setMethod("as.vector", signature(x = "diagonalMatrix", mode="missing"),
128              function(x, mode) {
129                  n <- x@Dim[1]
130                  mod <- mode(x@x)
131                  r <- vector(mod, length = n^2)
132                  if(n)
133                      r[1 + 0:(n - 1) * (n + 1)] <-
134                          if(x@diag == "U")
135                              switch(mod, "integer"= 1L,
136                                     "numeric"= 1, "logical"= TRUE)
137                          else x@x
138                  r
139              })
140    
141  setAs("diagonalMatrix", "generalMatrix", # prefer sparse:  setAs("diagonalMatrix", "generalMatrix", # prefer sparse:
142        function(from) as(from, paste(.M.kind(from), "gCMatrix", sep='')))        function(from) as(as(from, "CsparseMatrix"), "generalMatrix"))
143    
144  .diag.x <- function(m) {  .diag.x <- function(m) {
145      if(m@diag == "U")      if(m@diag == "U")
# Line 109  Line 153 
153      m      m
154  }  }
155    
156    if(FALSE) {
157  ## given the above, the following  4  coercions should be all unneeded;  ## given the above, the following  4  coercions should be all unneeded;
158  ## we prefer triangular to general:  ## we prefer triangular to general:
159  setAs("ddiMatrix", "dgTMatrix",  setAs("ddiMatrix", "dgTMatrix",
# Line 139  Line 184 
184    
185  setAs("ldiMatrix", "lgCMatrix",  setAs("ldiMatrix", "lgCMatrix",
186        function(from) as(as(from, "lgTMatrix"), "lgCMatrix"))        function(from) as(as(from, "lgTMatrix"), "lgCMatrix"))
187    }
188    
189    
190  if(FALSE) # now have faster  "ddense" -> "dge"  if(FALSE) # now have faster  "ddense" -> "dge"
# Line 189  Line 235 
235            function(x = 1, nrow, ncol) .diag.x(x))            function(x = 1, nrow, ncol) .diag.x(x))
236    
237    
238  subDiag <- function(x, i, j, drop) {  subDiag <- function(x, i, j, ..., drop) {
239      x <- as(x, "sparseMatrix")      x <- as(x, "sparseMatrix")
240      x <- if(missing(i))      x <- if(missing(i))
241          x[, j, drop=drop]          x[, j, drop=drop]
# Line 197  Line 243 
243          x[i, , drop=drop]          x[i, , drop=drop]
244      else      else
245          x[i,j, drop=drop]          x[i,j, drop=drop]
246      if(isDiagonal(x)) as(x, "diagonalMatrix") else x      if(isS4(x) && isDiagonal(x)) as(x, "diagonalMatrix") else x
247  }  }
248    
249  setMethod("[", signature(x = "diagonalMatrix", i = "index",  setMethod("[", signature(x = "diagonalMatrix", i = "index",
250                           j = "index", drop = "logical"), subDiag)                           j = "index", drop = "logical"), subDiag)
251  setMethod("[", signature(x = "diagonalMatrix", i = "index",  setMethod("[", signature(x = "diagonalMatrix", i = "index",
252                          j = "missing", drop = "logical"),                          j = "missing", drop = "logical"),
253            function(x, i, drop) subDiag(x, i=i, drop=drop))            function(x, i, j, ..., drop) subDiag(x, i=i, drop=drop))
254  setMethod("[", signature(x = "diagonalMatrix", i = "missing",  setMethod("[", signature(x = "diagonalMatrix", i = "missing",
255                           j = "index", drop = "logical"),                           j = "index", drop = "logical"),
256            function(x, j, drop) subDiag(x, j=j, drop=drop))            function(x, i, j, ..., drop) subDiag(x, j=j, drop=drop))
257    
258  ## When you assign to a diagonalMatrix, the result should be  ## When you assign to a diagonalMatrix, the result should be
259  ## diagonal or sparse ---  ## diagonal or sparse ---
260  ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch  ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch
261  replDiag <- function(x, i, j, value) {  ## Only(?) current bug:  x[i] <- value  is wrong when  i is *vector*
262    replDiag <- function(x, i, j, ..., value) {
263      x <- as(x, "sparseMatrix")      x <- as(x, "sparseMatrix")
     message(sprintf("diagnosing replDiag() -- nargs() = %d\n", nargs()))  
264      if(missing(i))      if(missing(i))
265          x[, j] <- value          x[, j] <- value
266      else if(missing(j))      else if(missing(j)) { ##  x[i , ] <- v  *OR*   x[i] <- v
267            na <- nargs()
268    ##         message("diagnosing replDiag() -- nargs()= ", na)
269            if(na == 4)
270          x[i, ] <- value          x[i, ] <- value
271      else          else if(na == 3)
272                x[i] <- value
273            else stop("Internal bug: nargs()=",na,"; please report")
274        } else
275          x[i,j] <- value          x[i,j] <- value
276      if(isDiagonal(x)) as(x, "diagonalMatrix") else x      if(isDiagonal(x)) as(x, "diagonalMatrix") else x
277  }  }
278    
279  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index",
280                                  j = "index", value = "replValue"), replDiag)                                  j = "index", value = "replValue"), replDiag)
281    
282  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index",
283                                  j = "missing", value = "replValue"),                                  j = "missing", value = "replValue"),
284                   function(x, i, value) replDiag(x, i=i, value=value))                   function(x,i,j, ..., value) {
285                         ## message("before replDiag() -- nargs()= ", nargs())
286                         if(nargs() == 3)
287                             replDiag(x, i=i, value=value)
288                         else ## nargs() == 4 :
289                             replDiag(x, i=i, , value=value)
290                     })
291    
292  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "matrix", # 2-col.matrix  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "matrix", # 2-col.matrix
293                                  j = "missing", value = "replValue"),                                  j = "missing", value = "replValue"),
294                   function(x, i, value) {                   function(x,i,j, ..., value) {
295                       if(ncol(i) == 2) {                       if(ncol(i) == 2) {
296                           if(all((ii <- i[,1]) == i[,2])) { # replace in diagonal only                           if(all((ii <- i[,1]) == i[,2])) { # replace in diagonal only
297                               x@x[ii] <- value                               x@x[ii] <- value
# Line 251  Line 311 
311    
312  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing",
313                                  j = "index", value = "replValue"),                                  j = "index", value = "replValue"),
314                   function(x, j, value) replDiag(x, j=j, value=value))                   function(x,i,j, ..., value) replDiag(x, j=j, value=value))
315    
316    
317  setMethod("t", signature(x = "diagonalMatrix"),  setMethod("t", signature(x = "diagonalMatrix"),
# Line 262  Line 322 
322  setMethod("isTriangular", signature(object = "diagonalMatrix"),  setMethod("isTriangular", signature(object = "diagonalMatrix"),
323            function(object) TRUE)            function(object) TRUE)
324  setMethod("isSymmetric", signature(object = "diagonalMatrix"),  setMethod("isSymmetric", signature(object = "diagonalMatrix"),
325            function(object) TRUE)            function(object, ...) TRUE)
326    
327    setMethod("symmpart", signature(x = "diagonalMatrix"), function(x) x)
328    setMethod("skewpart", signature(x = "diagonalMatrix"), setZero)
329    
330  setMethod("chol", signature(x = "ddiMatrix"),# pivot = "ANY"  setMethod("chol", signature(x = "ddiMatrix"),# pivot = "ANY"
331            function(x, pivot) {            function(x, pivot) {
# Line 366  Line 429 
429  ##           })  ##           })
430    
431  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
432            function(x, y = NULL) { x <- as(x, "sparseMatrix"); callGeneric() })            function(x, y = NULL) crossprod(as(x, "sparseMatrix"), y))
433    
434  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
435            function(x, y = NULL) { y <- as(y, "sparseMatrix"); callGeneric() })            function(x, y = NULL) crossprod(x, as(y, "sparseMatrix")))
436    
437  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
438            function(x, y = NULL) { x <- as(x, "sparseMatrix"); callGeneric() })            function(x, y = NULL) tcrossprod(as(x, "sparseMatrix"), y))
439    
440  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
441            function(x, y = NULL) { y <- as(y, "sparseMatrix"); callGeneric() })            function(x, y = NULL) tcrossprod(x, as(y, "sparseMatrix")))
442    
443    
444  ## 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 479 
479  setMethod("solve", signature(a = "diagonalMatrix", b = "Matrix"),  setMethod("solve", signature(a = "diagonalMatrix", b = "Matrix"),
480            solveDiag)            solveDiag)
481    
482    ## Schur()  ---> ./eigen.R
483    
484    
485    
# Line 456  Line 520 
520  ## For almost everything else, diag* shall be treated "as sparse" :  ## For almost everything else, diag* shall be treated "as sparse" :
521  ## These are cheap implementations via coercion  ## These are cheap implementations via coercion
522    
523  ## for disambiguation  ## for disambiguation --- define this for "sparseMatrix" , then for "ANY" :
524  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "sparseMatrix"),  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "sparseMatrix"),
525            function(e1,e2) callGeneric(as(e1, "sparseMatrix"), e2))            function(e1,e2) callGeneric(diag2tT.u(e1,e2), e2))
526  setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "diagonalMatrix"),  setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "diagonalMatrix"),
527            function(e1,e2) callGeneric(e1, as(e2, "sparseMatrix")))            function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1)))
528  ## in general:  ## in general:
529  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "ANY"),  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "ANY"),
530            function(e1,e2) callGeneric(as(e1,"sparseMatrix"), e2))            function(e1,e2) callGeneric(diag2tT.u(e1,e2), e2))
531  setMethod("Ops", signature(e1 = "ANY", e2 = "diagonalMatrix"),  setMethod("Ops", signature(e1 = "ANY", e2 = "diagonalMatrix"),
532            function(e1,e2) callGeneric(e1, as(e2,"sparseMatrix")))            function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1)))
533    
534    
535    

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

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