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 1805, Tue Mar 27 16:46:03 2007 UTC revision 2133, Sat Mar 15 18:56:50 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 43  Line 45 
45      ## make sure we had all matrices:      ## make sure we had all matrices:
46      if(!(is.matrix(dims) && nrow(dims) == 2))      if(!(is.matrix(dims) && nrow(dims) == 2))
47          stop("some arguments are not matrices")          stop("some arguments are not matrices")
48      csdim <- rbind(rep.int(0:0, 2),      csdim <- rbind(rep.int(0L, 2),
49                     apply(sapply(mlist, dim), 1, cumsum))                     apply(sapply(mlist, dim), 1, cumsum))
50      ret <- new("dgTMatrix", Dim = as.integer(csdim[nrow(csdim),]))      ret <- new("dgTMatrix", Dim = as.integer(csdim[nrow(csdim),]))
51      add1 <- matrix(1:0, 2,2)      add1 <- matrix(1:0, 2,2)
# Line 58  Line 60 
60      as(ret, "CsparseMatrix")      as(ret, "CsparseMatrix")
61  }  }
62    
63  diag2T <- function(from) {  
64      i <- if(from@diag == "U") integer(0) else seq_len(from@Dim[1]) - 1:1  .diag2tT <- function(from, uplo = "U") { ## to triangular Tsparse
65        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  setAs("diagonalMatrix", "triangularMatrix", diag2T)  diag2tT <- function(from) .diag2tT(from, "U")
74  setAs("diagonalMatrix", "sparseMatrix", diag2T)  
75    .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,
81            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)
98    setAs("diagonalMatrix", "sparseMatrix", diag2tT)
99  ## needed too (otherwise <dense> -> Tsparse is taken):  ## needed too (otherwise <dense> -> Tsparse is taken):
100  setAs("diagonalMatrix", "TsparseMatrix", diag2T)  setAs("diagonalMatrix", "TsparseMatrix", diag2tT)
101  ## is better than this:  ## is better than this:
102  ## setAs("diagonalMatrix", "sparseMatrix",  ## setAs("diagonalMatrix", "sparseMatrix",
103  ##       function(from)  ##       function(from)
104  ##        as(from, if(is(from, "dMatrix")) "dgCMatrix" else "lgCMatrix"))  ##        as(from, if(is(from, "dMatrix")) "dgCMatrix" else "lgCMatrix"))
105  setAs("diagonalMatrix", "CsparseMatrix",  setAs("diagonalMatrix", "CsparseMatrix",
106        function(from) as(diag2T(from), "CsparseMatrix"))        function(from) as(diag2tT(from), "CsparseMatrix"))
107    
108    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) {
# Line 85  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 100  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",
160        function(from) {        function(from) {
161            .Deprecated("as(, \"sparseMatrix\")")            .Deprecated("as(, \"sparseMatrix\")")
162            n <- from@Dim[1]            n <- from@Dim[1]
163            i <- seq_len(n) - 1:1            i <- seq_len(n) - 1L
164            new("dgTMatrix", i = i, j = i, x = .diag.x(from),            new("dgTMatrix", i = i, j = i, x = .diag.x(from),
165                Dim = c(n,n), Dimnames = from@Dimnames) })                Dim = c(n,n), Dimnames = from@Dimnames) })
166    
# Line 119  Line 173 
173            n <- from@Dim[1]            n <- from@Dim[1]
174            if(from@diag == "U") { # unit-diagonal            if(from@diag == "U") { # unit-diagonal
175                x <- rep.int(TRUE, n)                x <- rep.int(TRUE, n)
176                i <- seq_len(n) - 1:1                i <- seq_len(n) - 1L
177            } else { # "normal"            } else { # "normal"
178                nz <- nz.NA(from@x, na. = TRUE)                nz <- nz.NA(from@x, na. = TRUE)
179                x <- from@x[nz]                x <- from@x[nz]
180                i <- which(nz) - 1:1                i <- which(nz) - 1L
181            }            }
182            new("lgTMatrix", i = i, j = i, x = x,            new("lgTMatrix", i = i, j = i, x = x,
183                Dim = c(n,n), Dimnames = from@Dimnames) })                Dim = c(n,n), Dimnames = from@Dimnames) })
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 177  Line 232 
232    
233    
234  setMethod("diag", signature(x = "diagonalMatrix"),  setMethod("diag", signature(x = "diagonalMatrix"),
235            function(x = 1, nrow, ncol = n) .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 188  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")
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
293                                    j = "missing", value = "replValue"),
294                     function(x,i,j, ..., value) {
295                         if(ncol(i) == 2) {
296                             if(all((ii <- i[,1]) == i[,2])) { # replace in diagonal only
297                                 x@x[ii] <- value
298                                 x
299                             } else { ## no longer diagonal, but remain sparse:
300                                 x <- as(x, "sparseMatrix")
301                                 x[i] <- value
302                                 x
303                             }
304                         }
305                         else { # behave as "base R": use as if vector
306                             x <- as(x, "matrix")
307                             x[i] <- value
308                             Matrix(x)
309                         }
310                     })
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 232  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"),
331            function(x, pivot) {            function(x, pivot, ...) {
332                if(any(x@x < 0)) stop("chol() is undefined for diagonal matrix with negative entries")                if(any(x@x < 0))
333                      stop("chol() is undefined for diagonal matrix with negative entries")
334                x@x <- sqrt(x@x)                x@x <- sqrt(x@x)
335                x                x
336            })            })
337  ## chol(L) is L for logical diagonal:  ## chol(L) is L for logical diagonal:
338  setMethod("chol", signature(x = "ldiMatrix"), function(x, pivot) x)  setMethod("chol", signature(x = "ldiMatrix"), function(x, pivot, ...) x)
   
 setMethod("!", "ldiMatrix", function(e1) {  
     if(e1@diag == "N")  
         e1@x <- !e1@x  
     else { ## "U"  
         e1@diag <- "N"  
         e1@x <- rep.int(FALSE, e1@Dim[1])  
     }  
     e1  
 })  
339    
340  ## Basic Matrix Multiplication {many more to add}  ## Basic Matrix Multiplication {many more to add}
341  ##       ---------------------  ##       ---------------------
# Line 346  Line 430 
430  ##           })  ##           })
431    
432  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
433            function(x, y = NULL) { x <- as(x, "sparseMatrix"); callGeneric() })            function(x, y = NULL) crossprod(as(x, "sparseMatrix"), y))
434    
435  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
436            function(x, y = NULL) { y <- as(y, "sparseMatrix"); callGeneric() })            function(x, y = NULL) crossprod(x, as(y, "sparseMatrix")))
437    
438  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
439            function(x, y = NULL) { x <- as(x, "sparseMatrix"); callGeneric() })            function(x, y = NULL) tcrossprod(as(x, "sparseMatrix"), y))
440    
441  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
442            function(x, y = NULL) { y <- as(y, "sparseMatrix"); callGeneric() })            function(x, y = NULL) tcrossprod(x, as(y, "sparseMatrix")))
443    
444    
445  ## 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 396  Line 480 
480  setMethod("solve", signature(a = "diagonalMatrix", b = "Matrix"),  setMethod("solve", signature(a = "diagonalMatrix", b = "Matrix"),
481            solveDiag)            solveDiag)
482    
483    ## Schur()  ---> ./eigen.R
484    
485    
486    
# Line 428  Line 513 
513  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "ddiMatrix"),  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "ddiMatrix"),
514            diagOdiag)            diagOdiag)
515    
516    ## FIXME:    diagonal  o  triangular  |-->  triangular
517    ## -----     diagonal  o  symmetric   |-->  symmetric
518    ##    {also when other is sparse: do these "here" --
519    ##     before conversion to sparse, since that loses "diagonality"}
520    
521  ## For almost everything else, diag* shall be treated "as sparse" :  ## For almost everything else, diag* shall be treated "as sparse" :
522  ## These are cheap implementations via coercion  ## These are cheap implementations via coercion
523    
524  ## for disambiguation  ## for disambiguation --- define this for "sparseMatrix" , then for "ANY" :
525  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "sparseMatrix"),  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "sparseMatrix"),
526            function(e1,e2) callGeneric(as(e1, "sparseMatrix"), e2))            function(e1,e2) callGeneric(diag2tT.u(e1,e2), e2))
527  setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "diagonalMatrix"),  setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "diagonalMatrix"),
528            function(e1,e2) callGeneric(e1, as(e2, "sparseMatrix")))            function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1)))
529  ## in general:  ## in general:
530  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "ANY"),  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "ANY"),
531            function(e1,e2) callGeneric(as(e1,"sparseMatrix"), e2))            function(e1,e2) callGeneric(diag2tT.u(e1,e2), e2))
532  setMethod("Ops", signature(e1 = "ANY", e2 = "diagonalMatrix"),  setMethod("Ops", signature(e1 = "ANY", e2 = "diagonalMatrix"),
533            function(e1,e2) callGeneric(e1, as(e2,"sparseMatrix")))            function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1)))
534    
535    
536    

Legend:
Removed from v.1805  
changed lines
  Added in v.2133

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