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 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 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", kind = .M.kind(from)) {
65      new(paste(.M.kind(from), "tTMatrix", sep=''),      ## to triangular Tsparse
66        i <- if(from@diag == "U") integer(0) else seq_len(from@Dim[1]) - 1L
67        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  setAs("diagonalMatrix", "triangularMatrix", diag2T)  .diag2sT <- function(from, uplo = "U", kind = .M.kind(from)) {
75  setAs("diagonalMatrix", "sparseMatrix", diag2T)      ## to symmetric Tsparse
76        n <- from@Dim[1]
77        i <- seq_len(n) - 1L
78        new(paste(kind, "sTMatrix", sep=''),
79            Dim = from@Dim, Dimnames = from@Dimnames,
80            i = i, j = i, uplo = uplo,
81            x = if(from@diag == "N") from@x else ## "U"-diag
82            rep.int(switch(kind,
83                           "d" = 1.,
84                           "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", diag2T)  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(diag2T(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    
130    
131  setAs("diagonalMatrix", "matrix",  setAs("diagonalMatrix", "matrix",
132        function(from) {        function(from) {
# Line 85  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 100  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",
171        function(from) {        function(from) {
172            .Deprecated("as(, \"sparseMatrix\")")            .Deprecated("as(, \"sparseMatrix\")")
173            n <- from@Dim[1]            n <- from@Dim[1]
174            i <- seq_len(n) - 1:1            i <- seq_len(n) - 1L
175            new("dgTMatrix", i = i, j = i, x = .diag.x(from),            new("dgTMatrix", i = i, j = i, x = .diag.x(from),
176                Dim = c(n,n), Dimnames = from@Dimnames) })                Dim = c(n,n), Dimnames = from@Dimnames) })
177    
# Line 119  Line 184 
184            n <- from@Dim[1]            n <- from@Dim[1]
185            if(from@diag == "U") { # unit-diagonal            if(from@diag == "U") { # unit-diagonal
186                x <- rep.int(TRUE, n)                x <- rep.int(TRUE, n)
187                i <- seq_len(n) - 1:1                i <- seq_len(n) - 1L
188            } else { # "normal"            } else { # "normal"
189                nz <- nz.NA(from@x, na. = TRUE)                nz <- nz.NA(from@x, na. = TRUE)
190                x <- from@x[nz]                x <- from@x[nz]
191                i <- which(nz) - 1:1                i <- which(nz) - 1L
192            }            }
193            new("lgTMatrix", i = i, j = i, x = x,            new("lgTMatrix", i = i, j = i, x = x,
194                Dim = c(n,n), Dimnames = from@Dimnames) })                Dim = c(n,n), Dimnames = from@Dimnames) })
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 176  Line 242 
242        })        })
243    
244    
245  setMethod("diag", signature(x = "diagonalMatrix"),  ## In order to evade method dispatch ambiguity warnings,
246            function(x = 1, nrow, ncol = n) .diag.x(x))  ## 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))
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 188  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")
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
309                                    j = "missing", value = "replValue"),
310                     function(x,i,j, ..., value) {
311                         if(ncol(i) == 2) {
312                             if(all((ii <- i[,1]) == i[,2])) { # replace in diagonal only
313                                 x@x[ii] <- value
314                                 x
315                             } else { ## no longer diagonal, but remain sparse:
316                                 x <- as(x, "sparseMatrix")
317                                 x[i] <- value
318                                 x
319                             }
320                         }
321                         else { # behave as "base R": use as if vector
322                             x <- as(x, "matrix")
323                             x[i] <- value
324                             Matrix(x)
325                         }
326                     })
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 232  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)
   
 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  
 })  
355    
356  ## Basic Matrix Multiplication {many more to add}  ## Basic Matrix Multiplication {many more to add}
357  ##       ---------------------  ##       ---------------------
# Line 346  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 396  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 428  Line 529 
529  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "ddiMatrix"),  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "ddiMatrix"),
530            diagOdiag)            diagOdiag)
531    
532    ## FIXME:    diagonal  o  triangular  |-->  triangular
533    ## -----     diagonal  o  symmetric   |-->  symmetric
534    ##    {also when other is sparse: do these "here" --
535    ##     before conversion to sparse, since that loses "diagonality"}
536    
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.1805  
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