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 2197, Mon Jun 2 14:34:42 2008 UTC revision 2493, Mon Nov 2 11:55:05 2009 UTC
# Line 32  Line 32 
32      }      }
33  }  }
34    
35  ## Pkg 'spdep' had (relatively slow) versions of this as_dsCMatrix_I()  .sparseDiagonal <- function(n, x = rep.int(1,n), uplo = "U", shape = "t") {
 .symDiagonal <- function(n, x = rep.int(1,n), uplo = "U") {  
36      stopifnot(n == (n. <- as.integer(n)), (n <- n.) >= 0)      stopifnot(n == (n. <- as.integer(n)), (n <- n.) >= 0)
37      if((lx <- length(x)) == 1) x <- rep.int(x, n)      if((lx <- length(x)) == 1) x <- rep.int(x, n)
38      else if(lx != n) stop("length(x) must be 1 or n")      else if(lx != n) stop("length(x) must be 1 or n")
39      cls <-      stopifnot(is.character(shape), nchar(shape) == 1,
40          if(is.double(x)) "dsCMatrix"                any(shape == c("t","s","g"))) # triangular / symmetric / general
41          else if(is.logical(x)) "lsCMatrix"      kind <-
42            if(is.double(x)) "d"
43            else if(is.logical(x)) "l"
44          else { ## for now          else { ## for now
45              storage.mode(x) <- "double"              storage.mode(x) <- "double"
46              "dsCMatrix"              "d"
47          }          }
48      new(cls, Dim = c(n,n), x = x, uplo = uplo,      ii <- if(n) 0:(n - 1L) else integer(0)
49          i = if(n) 0:(n - 1L) else integer(0), p = 0:n)      if(shape == "g")
50            new(paste0(kind, "gCMatrix"), Dim = c(n,n),
51                x = x, i = ii, p = 0:n)
52        else new(paste0(kind, shape, "CMatrix"), Dim = c(n,n), uplo = uplo,
53                 x = x, i = ii, p = 0:n)
54  }  }
55    
56    ## Pkg 'spdep' had (relatively slow) versions of this as_dsCMatrix_I()
57    .symDiagonal <- function(n, x = rep.int(1,n), uplo = "U")
58        .sparseDiagonal(n, x, uplo, shape = "s")
59    
60    ## instead of   diagU2N(as(Diagonal(n), "CsparseMatrix")), diag = "N" in any case:
61    .trDiagonal <- function(n, x = rep.int(1,n), uplo = "U")
62        .sparseDiagonal(n, x, uplo, shape = "t")
63    
64    
65  ### 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.
66  ### 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
67  ### 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
68  ### who posted his bdiag() function written in December 1995.  ### who posted his bdiag() function written in December 1995.
69    if(FALSE)##--- no longer used:
70  bdiag <- function(...) {  .bdiag <- function(lst) {
71      if(nargs() == 0) return(new("dgCMatrix"))      ### block-diagonal matrix [a dgTMatrix] from list of matrices
72      ## else :      stopifnot(is.list(lst), length(lst) >= 1)
73      mlist <- if (nargs() == 1) as.list(...) else list(...)      dims <- sapply(lst, dim, USE.NAMES=FALSE)
     dims <- sapply(mlist, dim)  
74      ## make sure we had all matrices:      ## make sure we had all matrices:
75      if(!(is.matrix(dims) && nrow(dims) == 2))      if(!(is.matrix(dims) && nrow(dims) == 2))
76          stop("some arguments are not matrices")          stop("some arguments are not matrices")
77      csdim <- rbind(rep.int(0L, 2),      csdim <- rbind(rep.int(0L, 2),
78                     apply(sapply(mlist, dim), 1, cumsum))                     apply(dims, 1, cumsum))
79      ret <- new("dgTMatrix", Dim = as.integer(csdim[nrow(csdim),]))      r <- new("dgTMatrix")
80        r@Dim <- as.integer(csdim[nrow(csdim),])
81      add1 <- matrix(1:0, 2,2)      add1 <- matrix(1:0, 2,2)
82      for(i in seq_along(mlist)) {      for(i in seq_along(lst)) {
83          indx <- apply(csdim[i:(i+1),] + add1, 2, function(n) n[1]:n[2])          indx <- apply(csdim[i:(i+1),] + add1, 2, function(n) n[1]:n[2])
84          if(is.null(dim(indx))) ## non-square matrix          if(is.null(dim(indx))) ## non-square matrix
85              ret[indx[[1]],indx[[2]]] <- mlist[[i]]              r[indx[[1]],indx[[2]]] <- lst[[i]]
86          else ## square matrix          else ## square matrix
87              ret[indx[,1],indx[,2]] <- mlist[[i]]              r[indx[,1], indx[,2]] <- lst[[i]]
88        }
89        r
90    }
91    ## expand(<mer>) needed something like bdiag() for lower-triangular
92    ## (Tsparse) Matrices; hence Doug Bates provided a much more efficient
93    ##  implementation for those; now extended and generalized:
94    .bdiag <- function(lst) {
95        ## block-diagonal matrix [a dgTMatrix] from list of matrices
96        stopifnot(is.list(lst), (nl <- length(lst)) >= 1)
97    
98        Tlst <- lapply(lapply(lst, Matrix:::as_Csp2), # includes "diagU2N"
99                       as, "TsparseMatrix")
100        if(nl == 1) return(Tlst[[1]])
101        ## else
102        i_off <- c(0L, cumsum(sapply(Tlst, nrow)))
103        j_off <- c(0L, cumsum(sapply(Tlst, ncol)))
104    
105        clss <- sapply(Tlst, class)
106        knds <- substr(clss, 2, 2)
107        sym  <- knds == "s" # symmetric ones
108        tri  <- knds == "t" # triangular ones
109        use.n <- any(is.n <- substr(clss,1,1) == "n")
110        if(use.n && !(use.n <- all(is.n)))
111            Tlst[is.n] <- lapply(Tlst[is.n], as, "lMatrix")
112        if(all(sym)) { ## result should be *symmetric*
113            uplos <- sapply(Tlst, slot, "uplo") ## either "U" or "L"
114            tLU <- table(uplos)# of length 1 or 2 ..
115            if(length(tLU) == 1) { ## all "U" or all "L"
116                useU <- uplos[1] == "U"
117            } else { ## length(tLU) == 2, counting "L" and "U"
118                useU <- diff(tLU) >= 0
119                if(useU && (hasL <- tLU[1] > 0))
120                    Tlst[hasL] <- lapply(Tlst[hasL], t)
121                else if(!useU && (hasU <- tLU[2] > 0))
122                    Tlst[hasU] <- lapply(Tlst[hasU], t)
123            }
124            if(use.n) { ## return nsparseMatrix :
125                r <- new("nsTMatrix")
126            } else {
127                r <- new("dsTMatrix")
128                r@x <- unlist(lapply(Tlst, slot, "x"))
129            }
130            r@uplo <- if(useU) "U" else "L"
131        }
132        else if(all(tri) && { ULs <- sapply(Tlst, slot, "uplo")##  "U" or "L"
133                              all(ULs[1L] == ULs[-1L]) } ## all upper or all lower
134           ){ ## *triangular* result
135    
136            if(use.n) { ## return nsparseMatrix :
137                r <- new("ntTMatrix")
138            } else {
139                r <- new("dtTMatrix")
140                r@x <- unlist(lapply(Tlst, slot, "x"))
141            }
142            r@uplo <- ULs[1L]
143        }
144        else {
145            if(any(sym))
146                Tlst[sym] <- lapply(Tlst[sym], as, "generalMatrix")
147            if(use.n) { ## return nsparseMatrix :
148                r <- new("ngTMatrix")
149            } else {
150                r <- new("dgTMatrix")
151                r@x <- unlist(lapply(Tlst, slot, "x"))
152            }
153      }      }
154      ## slightly debatable if we really should return Csparse.. :      r@Dim <- c(i_off[nl+1], j_off[nl + 1])
155      as(ret, "CsparseMatrix")      r@i <- unlist(lapply(1:nl, function(k) Tlst[[k]]@i + i_off[k]))
156        r@j <- unlist(lapply(1:nl, function(k) Tlst[[k]]@j + j_off[k]))
157        r
158    }
159    
160    bdiag <- function(...) {
161        if((nA <- nargs()) == 0) return(new("dgCMatrix"))
162        if(nA == 1 && !is.list(...))
163            return(as(..., "CsparseMatrix"))
164        alis <- if(nA == 1 && is.list(..1)) ..1 else list(...)
165        if(length(alis) == 1)
166            return(as(alis[[1]], "CsparseMatrix"))
167    
168        ## else : two or more arguments
169        as(.bdiag(alis), "CsparseMatrix")
170  }  }
171    
172    
# Line 124  Line 217 
217  ## ddi*:  ## ddi*:
218  diag2tT <- function(from) .diag2tT(from, "U", "d")  diag2tT <- function(from) .diag2tT(from, "U", "d")
219  setAs("ddiMatrix", "triangularMatrix", diag2tT)  setAs("ddiMatrix", "triangularMatrix", diag2tT)
220  setAs("ddiMatrix", "sparseMatrix", diag2tT)  ##_no_longer_ setAs("ddiMatrix", "sparseMatrix", diag2tT)
221  ## needed too (otherwise <dense> -> Tsparse is taken):  ## needed too (otherwise <dense> -> Tsparse is taken):
222  setAs("ddiMatrix", "TsparseMatrix", diag2tT)  setAs("ddiMatrix", "TsparseMatrix", diag2tT)
223    setAs("ddiMatrix", "dsparseMatrix", diag2tT)
224  setAs("ddiMatrix", "CsparseMatrix",  setAs("ddiMatrix", "CsparseMatrix",
225        function(from) as(.diag2tT(from, "U", "d"), "CsparseMatrix"))        function(from) as(.diag2tT(from, "U", "d"), "CsparseMatrix"))
226  setAs("ddiMatrix", "symmetricMatrix",  setAs("ddiMatrix", "symmetricMatrix",
# Line 135  Line 229 
229  ## ldi*:  ## ldi*:
230  diag2tT <- function(from) .diag2tT(from, "U", "l")  diag2tT <- function(from) .diag2tT(from, "U", "l")
231  setAs("ldiMatrix", "triangularMatrix", diag2tT)  setAs("ldiMatrix", "triangularMatrix", diag2tT)
232  setAs("ldiMatrix", "sparseMatrix", diag2tT)  ##_no_longer_ setAs("ldiMatrix", "sparseMatrix", diag2tT)
233  ## needed too (otherwise <dense> -> Tsparse is taken):  ## needed too (otherwise <dense> -> Tsparse is taken):
234  setAs("ldiMatrix", "TsparseMatrix", diag2tT)  setAs("ldiMatrix", "TsparseMatrix", diag2tT)
235    setAs("ldiMatrix", "lsparseMatrix", diag2tT)
236  setAs("ldiMatrix", "CsparseMatrix",  setAs("ldiMatrix", "CsparseMatrix",
237        function(from) as(.diag2tT(from, "U", "l"), "CsparseMatrix"))        function(from) as(.diag2tT(from, "U", "l"), "CsparseMatrix"))
238  setAs("ldiMatrix", "symmetricMatrix",  setAs("ldiMatrix", "symmetricMatrix",
# Line 152  Line 247 
247                Dim = from@Dim, Dimnames = from@Dimnames)                Dim = from@Dim, Dimnames = from@Dimnames)
248        })        })
249    
250    setAs("diagonalMatrix", "nsparseMatrix", function(from) as(from, "nMatrix"))
251    
252  ## Cheap fast substitute for diag() which *does* preserve the mode of x :  ## Cheap fast substitute for diag() which *does* preserve the mode of x :
253  mkDiag <- function(x, n) {  mkDiag <- function(x, n) {
# Line 181  Line 277 
277  setAs("diagonalMatrix", "generalMatrix", # prefer sparse:  setAs("diagonalMatrix", "generalMatrix", # prefer sparse:
278        function(from) as(as(from, "CsparseMatrix"), "generalMatrix"))        function(from) as(as(from, "CsparseMatrix"), "generalMatrix"))
279    
280    setAs("diagonalMatrix", "denseMatrix",
281          function(from) as(as(from, "CsparseMatrix"), "denseMatrix"))
282    
283  .diag.x <- function(m) if(m@diag == "U") rep.int(as1(m@x), m@Dim[1]) else m@x  .diag.x <- function(m) if(m@diag == "U") rep.int(as1(m@x), m@Dim[1]) else m@x
284    
285  .diag.2N <- function(m) {  .diag.2N <- function(m) {
# Line 201  Line 300 
300            d <- dim(from)            d <- dim(from)
301            if(d[1] != (n <- d[2])) stop("non-square matrix")            if(d[1] != (n <- d[2])) stop("non-square matrix")
302            if(any(from[row(from) != col(from)] != 0))            if(any(from[row(from) != col(from)] != 0))
303                stop("has non-zero off-diagonal entries")                stop("matrix with non-zero off-diagonals cannot be coerced to diagonalMatrix")
304            x <- diag(from)            x <- diag(from)
305            if(is.logical(x)) {            if(is.logical(x)) {
306                cl <- "ldiMatrix"                cl <- "ldiMatrix"
307                uni <- all(x)                uni <- allTrue(x) ## uni := {is it unit-diagonal ?}
308            } else {            } else {
309                cl <- "ddiMatrix"                cl <- "ddiMatrix"
310                uni <- all(x == 1)                uni <- allTrue(x == 1)
311                storage.mode(x) <- "double"                storage.mode(x) <- "double"
312            } ## TODO: complex            } ## TODO: complex
313            new(cl, Dim = c(n,n), diag = if(uni) "U" else "N",            new(cl, Dim = c(n,n), diag = if(uni) "U" else "N",
# Line 225  Line 324 
324            x <- diag(from)            x <- diag(from)
325            if(is.logical(x)) {            if(is.logical(x)) {
326                cl <- "ldiMatrix"                cl <- "ldiMatrix"
327                uni <- all(x)                uni <- allTrue(x)
328            } else {            } else {
329                cl <- "ddiMatrix"                cl <- "ddiMatrix"
330                uni <- all(x == 1)                uni <- allTrue(x == 1)
331                storage.mode(x) <- "double"                storage.mode(x) <- "double"
332            }            } ## TODO: complex
333            new(cl, Dim = c(n,n), diag = if(uni) "U" else "N",            new(cl, Dim = c(n,n), diag = if(uni) "U" else "N",
334                x = if(uni) x[FALSE] else x)                x = if(uni) x[FALSE] else x)
335        })        })
# Line 240  Line 339 
339            function(x = 1, nrow, ncol) .diag.x(x))            function(x = 1, nrow, ncol) .diag.x(x))
340    
341  subDiag <- function(x, i, j, ..., drop) {  subDiag <- function(x, i, j, ..., drop) {
342      x <- as(x, "sparseMatrix")      x <- as(x, "TsparseMatrix")
343      x <- if(missing(i))      x <- if(missing(i))
344          x[, j, drop=drop]          x[, j, drop=drop]
345      else if(missing(j))      else if(missing(j))
346          x[i, , drop=drop]          if(nargs() == 4) x[i, , drop=drop] else x[i, drop=drop]
347      else      else
348          x[i,j, drop=drop]          x[i,j, drop=drop]
349      if(isS4(x) && isDiagonal(x)) as(x, "diagonalMatrix") else x      if(isS4(x) && isDiagonal(x)) as(x, "diagonalMatrix") else x
# Line 254  Line 353 
353                           j = "index", drop = "logical"), subDiag)                           j = "index", drop = "logical"), subDiag)
354  setMethod("[", signature(x = "diagonalMatrix", i = "index",  setMethod("[", signature(x = "diagonalMatrix", i = "index",
355                          j = "missing", drop = "logical"),                          j = "missing", drop = "logical"),
356            function(x, i, j, ..., drop) subDiag(x, i=i, drop=drop))            function(x, i, j, ..., drop) {
357                  na <- nargs()
358                  Matrix.msg("diag[i,m,l] : nargs()=", na, .M.level = 2)
359                  if(na == 4)
360                       subDiag(x, i=i, , drop=drop)
361                  else subDiag(x, i=i,   drop=drop)
362              })
363  setMethod("[", signature(x = "diagonalMatrix", i = "missing",  setMethod("[", signature(x = "diagonalMatrix", i = "missing",
364                           j = "index", drop = "logical"),                           j = "index", drop = "logical"),
365            function(x, i, j, ..., drop) subDiag(x, j=j, drop=drop))            function(x, i, j, ..., drop) subDiag(x, j=j, drop=drop))
# Line 264  Line 369 
369  ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch  ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch
370  ## Only(?) current bug:  x[i] <- value  is wrong when  i is *vector*  ## Only(?) current bug:  x[i] <- value  is wrong when  i is *vector*
371  replDiag <- function(x, i, j, ..., value) {  replDiag <- function(x, i, j, ..., value) {
372      x <- as(x, "sparseMatrix")      x <- as(x, "TsparseMatrix")
373      if(missing(i))      if(missing(i))
374          x[, j] <- value          x[, j] <- value
375      else if(missing(j)) { ##  x[i , ] <- v  *OR*   x[i] <- v      else if(missing(j)) { ##  x[i , ] <- v  *OR*   x[i] <- v
# Line 309  Line 414 
414                               x@x[ii] <- value                               x@x[ii] <- value
415                               x                               x
416                           } else { ## no longer diagonal, but remain sparse:                           } else { ## no longer diagonal, but remain sparse:
417                               x <- as(x, "sparseMatrix")                               x <- as(x, "TsparseMatrix")
418                               x[i] <- value                               x[i] <- value
419                               x                               x
420                           }                           }
# Line 325  Line 430 
430                                  j = "index", value = "replValue"),                                  j = "index", value = "replValue"),
431                   function(x,i,j, ..., value) replDiag(x, j=j, value=value))                   function(x,i,j, ..., value) replDiag(x, j=j, value=value))
432    
433    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",
434                                    value = "sparseMatrix"),
435                     function (x, i, j, ..., value)
436                     callGeneric(x=x, , j=j, value = as(value, "sparseVector")))
437    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "missing",
438                                    value = "sparseMatrix"),
439                     function (x, i, j, ..., value)
440                     callGeneric(x=x, i=i, , value = as(value, "sparseVector")))
441    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "index",
442                                    value = "sparseMatrix"),
443                     function (x, i, j, ..., value)
444                     callGeneric(x=x, i=i, j=j, value = as(value, "sparseVector")))
445    
446    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",
447                                    value = "sparseVector"),
448                     replDiag)
449    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "missing",
450                                    value = "sparseVector"),
451                     replDiag)
452    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "index",
453                                    value = "sparseVector"),
454                     replDiag)
455    
456    
457  setMethod("t", signature(x = "diagonalMatrix"),  setMethod("t", signature(x = "diagonalMatrix"),
458            function(x) { x@Dimnames <- x@Dimnames[2:1] ; x })            function(x) { x@Dimnames <- x@Dimnames[2:1] ; x })
# Line 472  Line 600 
600  ##        function(x, y = NULL) {  ##        function(x, y = NULL) {
601  ##           })  ##           })
602    
603  ## setMethod("tcrossprod", signature(x = "denseMatrix", y = "diagonalMatrix"),  Cspdiagprod <- function(x, y) {
604  ##        function(x, y = NULL) {      dx <- dim(x <- .Call(Csparse_diagU2N, x))
605  ##           })      dy <- dim(y)
606        if(dx[2] != dy[1]) stop("non-matching dimensions")
607        ind <- rep.int(seq_len(dx[2]), x@p[-1] - x@p[-dx[2]-1L])
608        if(y@diag == "N")
609            x@x <- x@x * y@x[ind]
610        x
611    }
612    
613    diagCspprod <- function(x, y) {
614        dx <- dim(x)
615        dy <- dim(y <- .Call(Csparse_diagU2N, y))
616        if(dx[2] != dy[1]) stop("non-matching dimensions")
617        if(x@diag == "N")
618            y@x <- y@x * x@x[y@i + 1L]
619        y
620    }
621    
622    setMethod("crossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
623              function(x, y = NULL) diagCspprod(x, y))
624    
625  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
626            function(x, y = NULL) crossprod(as(x, "sparseMatrix"), y))            function(x, y = NULL) diagCspprod(x, as(y, "CsparseMatrix")))
627    
628    ## Prefer calling diagCspprod to Cspdiagprod if going to transpose anyway
629    ##  x'y == (y'x)'
630    setMethod("crossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
631              function(x, y = NULL) t(diagCspprod(y, x)))
632    
633  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
634            function(x, y = NULL) crossprod(x, as(y, "sparseMatrix")))            function(x, y = NULL) t(diagCspprod(y, as(x, "Csparsematrix"))))
635    
636    setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
637              function(x, y = NULL) diagCspprod(x, t(y)))
638    
639  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
640            function(x, y = NULL) tcrossprod(as(x, "sparseMatrix"), y))            function(x, y = NULL) diagCspprod(x, t(as(y, "CsparseMatrix"))))
641    
642    setMethod("tcrossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
643              function(x, y = NULL) Cspdiagprod(x, y))
644    
645  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
646            function(x, y = NULL) tcrossprod(x, as(y, "sparseMatrix")))            function(x, y = NULL) Cspdiagprod(as(x, "CsparseMatrix"), y))
647    
648    setMethod("%*%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
649              function(x, y) diagCspprod(x, y))
650    
 ## FIXME?: In theory, this can be done *FASTER*, in some cases, via tapply1()  
651  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
652            function(x, y) as(x, "sparseMatrix") %*% y)            function(x, y) diagCspprod(as(x, "CsparseMatrix"), y))
653    
654  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
655            function(x, y) x %*% as(y, "sparseMatrix"))            function(x, y) Cspdiagprod(as(x, "CsparseMatrix"), y))
656  ## NB: The previous is *not* triggering for  "ddi" o "dgC" (= distance 3)  
 ##     since there's a "ddense" o "Csparse" at dist. 2 => triggers first.  
 ## ==> do this:  
 setMethod("%*%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),  
           function(x, y) as(x, "CsparseMatrix") %*% y)  
657  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
658            function(x, y) x %*% as(y, "CsparseMatrix"))            function(x, y) Cspdiagprod(x, y))
659  ## NB: this is *not* needed for Tsparse & Rsparse  
660  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*
661  ##       do indeed work by going through sparse (and *not* ddense)!  ##       do indeed work by going through sparse (and *not* ddense)!
662    
   
   
663  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),
664            function(a, b, ...) {            function(a, b, ...) {
665                a@x <- 1/ a@x                a@x <- 1/ a@x
# Line 559  Line 712 
712          isNum <- (is.numeric(r) || is.numeric(r00))          isNum <- (is.numeric(r) || is.numeric(r00))
713          isLog <- (is.logical(r) || is.logical(r00))          isLog <- (is.logical(r) || is.logical(r00))
714    
715          if(getOption("verbose"))          Matrix.msg("exploding   <diag>  o  <diag>  into dense matrix")
             message("exploding  <diag>  o  <diag>  into dense matrix")  
716          d <- e1@Dim          d <- e1@Dim
717          n <- d[1]          n <- d[1]
718          stopifnot(length(r) == n)          stopifnot(length(r) == n)
# Line 588  Line 740 
740              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)
741  }  }
742    
743  ## FIXME:    diagonal  o  triangular  |-->  triangular  ## diagonal  o  triangular  |-->  triangular
744  ## -----     diagonal  o  symmetric   |-->  symmetric  ## diagonal  o  symmetric   |-->  symmetric
745  ##    {also when other is sparse: do these "here" --  ##    {also when other is sparse: do these "here" --
746  ##     before conversion to sparse, since that loses "diagonality"}  ##     before conversion to sparse, since that loses "diagonality"}
747    
# Line 619  Line 771 
771  ## ddi*: Arith: result numeric, potentially ddiMatrix  ## ddi*: Arith: result numeric, potentially ddiMatrix
772  setMethod("Arith", signature(e1 = "ddiMatrix", e2 = "numeric"),  setMethod("Arith", signature(e1 = "ddiMatrix", e2 = "numeric"),
773            function(e1,e2) {            function(e1,e2) {
774                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]; nsq <- n^2
775                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
776                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
777                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
# Line 637  Line 789 
789    
790  setMethod("Arith", signature(e1 = "numeric", e2 = "ddiMatrix"),  setMethod("Arith", signature(e1 = "numeric", e2 = "ddiMatrix"),
791            function(e1,e2) {            function(e1,e2) {
792                n <- e2@Dim[1]; nsq <- n*n                n <- e2@Dim[1]; nsq <- n^2
793                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
794                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
795                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
# Line 656  Line 808 
808  ## ldi* Arith --> result numeric, potentially ddiMatrix  ## ldi* Arith --> result numeric, potentially ddiMatrix
809  setMethod("Arith", signature(e1 = "ldiMatrix", e2 = "numeric"),  setMethod("Arith", signature(e1 = "ldiMatrix", e2 = "numeric"),
810            function(e1,e2) {            function(e1,e2) {
811                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]; nsq <- n^2
812                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
813                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
814                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
# Line 675  Line 827 
827    
828  setMethod("Arith", signature(e1 = "numeric", e2 = "ldiMatrix"),  setMethod("Arith", signature(e1 = "numeric", e2 = "ldiMatrix"),
829            function(e1,e2) {            function(e1,e2) {
830                n <- e2@Dim[1]; nsq <- n*n                n <- e2@Dim[1]; nsq <- n^2
831                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
832                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
833                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
# Line 695  Line 847 
847  ## ddi*: for "Ops" without Arith --> result logical, potentially ldi  ## ddi*: for "Ops" without Arith --> result logical, potentially ldi
848  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "numeric"),  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "numeric"),
849            function(e1,e2) {            function(e1,e2) {
850                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]; nsq <- n^2
851                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
852                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
853                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
# Line 714  Line 866 
866    
867  setMethod("Ops", signature(e1 = "numeric", e2 = "ddiMatrix"),  setMethod("Ops", signature(e1 = "numeric", e2 = "ddiMatrix"),
868            function(e1,e2) {            function(e1,e2) {
869                n <- e2@Dim[1]; nsq <- n*n                n <- e2@Dim[1]; nsq <- n^2
870                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
871                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
872                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
# Line 734  Line 886 
886  ## ldi*: for "Ops" without Arith --> result logical, potentially ldi  ## ldi*: for "Ops" without Arith --> result logical, potentially ldi
887  setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "numeric"),  setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "numeric"),
888            function(e1,e2) {            function(e1,e2) {
889                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]; nsq <- n^2
890                f0 <- callGeneric(FALSE, e2)                f0 <- callGeneric(FALSE, e2)
891                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
892                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
# Line 752  Line 904 
904    
905  setMethod("Ops", signature(e1 = "numeric", e2 = "ldiMatrix"),  setMethod("Ops", signature(e1 = "numeric", e2 = "ldiMatrix"),
906            function(e1,e2) {            function(e1,e2) {
907                n <- e2@Dim[1]; nsq <- n*n                n <- e2@Dim[1]; nsq <- n^2
908                f0 <- callGeneric(e1, FALSE)                f0 <- callGeneric(e1, FALSE)
909                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
910                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
# Line 774  Line 926 
926  for(other in c("ANY", "Matrix", "dMatrix")) {  for(other in c("ANY", "Matrix", "dMatrix")) {
927      ## ddi*:      ## ddi*:
928      setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other),      setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other),
929                function(e1,e2) callGeneric(diag2tT.u(e1,e2, "d"), e2))                function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2))
930      setMethod("Ops", signature(e1 = other, e2 = "ddiMatrix"),      setMethod("Ops", signature(e1 = other, e2 = "ddiMatrix"),
931                function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "d")))                function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "d")))
932      ## ldi*:      ## ldi*:
933      setMethod("Ops", signature(e1 = "ldiMatrix", e2 = other),      setMethod("Ops", signature(e1 = "ldiMatrix", e2 = other),
934                function(e1,e2) callGeneric(diag2tT.u(e1,e2, "l"), e2))                function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2))
935      setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),      setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),
936                function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "l")))                function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "l")))
937    }
938    
939    ## Direct subclasses of "denseMatrix": currently ddenseMatrix, ldense... :
940    dense.subCl <- local({ dM.scl <- getClass("denseMatrix")@subclasses
941                           names(dM.scl)[sapply(dM.scl, slot, "distance") == 1] })
942    for(DI in diCls) {
943        for(c2 in c(dense.subCl, "Matrix")) {
944            for(Fun in c("*", "^", "&")) {
945                setMethod(Fun, signature(e1 = DI, e2 = c2),
946                          function(e1,e2) callGeneric(e1, Diagonal(x = diag(e2))))
947                setMethod(Fun, signature(e1 = c2, e2 = DI),
948                          function(e1,e2) callGeneric(Diagonal(x = diag(e1)), e2))
949            }
950            ## NB: This arguably implicitly uses  0/0 :== 0  to keep diagonality
951            for(Fun in c("%%", "%/%", "/")) {
952                setMethod(Fun, signature(e1 = DI, e2 = c2),
953                          function(e1,e2) callGeneric(e1, Diagonal(x = diag(e2))))
954  }  }
 ## This should *not* dispatch to <dense> methods (in ./Ops.R ), as  
 ##  FALSE & <anything> |-> FALSE : hence result should be diagonal:  
 for(cl in diCls) {  
     setMethod("&", signature(e1 = cl, e2 = "ANY"),  
               function(e1,e2) e1 & as(e2,"Matrix"))  
     setMethod("&", signature(e1 = "ANY", e2 = cl),  
               function(e1,e2) as(e1,"Matrix") & e2)  
     for(c2 in c("denseMatrix", "Matrix")) {  
         setMethod("&", signature(e1 = cl, e2 = c2),  
                   function(e1,e2) e1 & Diagonal(x = diag(e2)))  
         setMethod("&", signature(e1 = c2, e2 = cl),  
                   function(e1,e2) Diagonal(x = diag(e1)) & e2)  
955      }      }
956  }  }
957    
958    
959  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"
960  ### ----------  any, all: separately here  ### ----------   the last 4: separately here
961  for(cl in diCls) {  for(cl in diCls) {
962  setMethod("any", cl,  setMethod("any", cl,
963            function (x, ..., na.rm) {            function (x, ..., na.rm) {
964                if(any(x@Dim == 0)) FALSE                if(any(x@Dim == 0)) FALSE
965                else if(x@diag == "U") TRUE else any(x@x, ..., na.rm = na.rm)                else if(x@diag == "U") TRUE else any(x@x, ..., na.rm = na.rm)
966            })            })
967  setMethod("all",  cl, function (x, ..., na.rm) any(x@Dim == 0))  setMethod("all",  cl, function (x, ..., na.rm) {
968  setMethod("prod", cl, function (x, ..., na.rm) as.numeric(any(x@Dim == 0)))      n <- x@Dim[1]
969        if(n >= 2) FALSE
970        else if(n == 0 || x@diag == "U") TRUE
971        else all(x@x, ..., na.rm = na.rm)
972    })
973    setMethod("prod", cl, function (x, ..., na.rm) {
974        n <- x@Dim[1]
975        if(n >= 2) 0
976        else if(n == 0 || x@diag == "U") 1
977        else ## n == 1, diag = "N" :
978            prod(x@x, ..., na.rm = na.rm)
979    })
980    
981  setMethod("sum", cl,  setMethod("sum", cl,
982            function(x, ..., na.rm) {            function(x, ..., na.rm) {
# Line 847  Line 1015 
1015      invisible(x)      invisible(x)
1016  }  }
1017    
1018    ## somewhat consistent with "print" for sparseMatrix :
1019    setMethod("print", signature(x = "diagonalMatrix"), prDiag)
1020    
1021  setMethod("show", signature(object = "diagonalMatrix"),  setMethod("show", signature(object = "diagonalMatrix"),
1022            function(object) {            function(object) {
1023                d <- dim(object)                d <- dim(object)
1024                cl <- class(object)                cl <- class(object)
1025                cat(sprintf('%d x %d diagonal matrix of class "%s"\n',                cat(sprintf('%d x %d diagonal matrix of class "%s"',
1026                            d[1], d[2], cl))                            d[1], d[2], cl))
1027                  if(d[1] < 50) {
1028                      cat("\n")
1029                prDiag(object)                prDiag(object)
1030                  } else {
1031                      cat(", with diagonal entries\n")
1032                      show(diag(object))
1033                      invisible(object)
1034                  }
1035            })            })

Legend:
Removed from v.2197  
changed lines
  Added in v.2493

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