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 2472, Sat Sep 19 06:10:48 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      ## slightly debatable if we really should return Csparse.. :      }
154      as(ret, "CsparseMatrix")      r@Dim <- c(i_off[nl+1], j_off[nl + 1])
155        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", "CsparseMatrix",  setAs("ddiMatrix", "CsparseMatrix",
# Line 135  Line 228 
228  ## ldi*:  ## ldi*:
229  diag2tT <- function(from) .diag2tT(from, "U", "l")  diag2tT <- function(from) .diag2tT(from, "U", "l")
230  setAs("ldiMatrix", "triangularMatrix", diag2tT)  setAs("ldiMatrix", "triangularMatrix", diag2tT)
231  setAs("ldiMatrix", "sparseMatrix", diag2tT)  ##_no_longer_ setAs("ldiMatrix", "sparseMatrix", diag2tT)
232  ## needed too (otherwise <dense> -> Tsparse is taken):  ## needed too (otherwise <dense> -> Tsparse is taken):
233  setAs("ldiMatrix", "TsparseMatrix", diag2tT)  setAs("ldiMatrix", "TsparseMatrix", diag2tT)
234  setAs("ldiMatrix", "CsparseMatrix",  setAs("ldiMatrix", "CsparseMatrix",
# Line 152  Line 245 
245                Dim = from@Dim, Dimnames = from@Dimnames)                Dim = from@Dim, Dimnames = from@Dimnames)
246        })        })
247    
248    setAs("diagonalMatrix", "nsparseMatrix", function(from) as(from, "nMatrix"))
249    
250  ## Cheap fast substitute for diag() which *does* preserve the mode of x :  ## Cheap fast substitute for diag() which *does* preserve the mode of x :
251  mkDiag <- function(x, n) {  mkDiag <- function(x, n) {
# Line 181  Line 275 
275  setAs("diagonalMatrix", "generalMatrix", # prefer sparse:  setAs("diagonalMatrix", "generalMatrix", # prefer sparse:
276        function(from) as(as(from, "CsparseMatrix"), "generalMatrix"))        function(from) as(as(from, "CsparseMatrix"), "generalMatrix"))
277    
278    setAs("diagonalMatrix", "denseMatrix",
279          function(from) as(as(from, "CsparseMatrix"), "denseMatrix"))
280    
281  .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
282    
283  .diag.2N <- function(m) {  .diag.2N <- function(m) {
# Line 201  Line 298 
298            d <- dim(from)            d <- dim(from)
299            if(d[1] != (n <- d[2])) stop("non-square matrix")            if(d[1] != (n <- d[2])) stop("non-square matrix")
300            if(any(from[row(from) != col(from)] != 0))            if(any(from[row(from) != col(from)] != 0))
301                stop("has non-zero off-diagonal entries")                stop("matrix with non-zero off-diagonals cannot be coerced to diagonalMatrix")
302            x <- diag(from)            x <- diag(from)
303            if(is.logical(x)) {            if(is.logical(x)) {
304                cl <- "ldiMatrix"                cl <- "ldiMatrix"
305                uni <- all(x)                uni <- allTrue(x) ## uni := {is it unit-diagonal ?}
306            } else {            } else {
307                cl <- "ddiMatrix"                cl <- "ddiMatrix"
308                uni <- all(x == 1)                uni <- allTrue(x == 1)
309                storage.mode(x) <- "double"                storage.mode(x) <- "double"
310            } ## TODO: complex            } ## TODO: complex
311            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 322 
322            x <- diag(from)            x <- diag(from)
323            if(is.logical(x)) {            if(is.logical(x)) {
324                cl <- "ldiMatrix"                cl <- "ldiMatrix"
325                uni <- all(x)                uni <- allTrue(x)
326            } else {            } else {
327                cl <- "ddiMatrix"                cl <- "ddiMatrix"
328                uni <- all(x == 1)                uni <- allTrue(x == 1)
329                storage.mode(x) <- "double"                storage.mode(x) <- "double"
330            }            } ## TODO: complex
331            new(cl, Dim = c(n,n), diag = if(uni) "U" else "N",            new(cl, Dim = c(n,n), diag = if(uni) "U" else "N",
332                x = if(uni) x[FALSE] else x)                x = if(uni) x[FALSE] else x)
333        })        })
# Line 240  Line 337 
337            function(x = 1, nrow, ncol) .diag.x(x))            function(x = 1, nrow, ncol) .diag.x(x))
338    
339  subDiag <- function(x, i, j, ..., drop) {  subDiag <- function(x, i, j, ..., drop) {
340      x <- as(x, "sparseMatrix")      x <- as(x, "TsparseMatrix")
341      x <- if(missing(i))      x <- if(missing(i))
342          x[, j, drop=drop]          x[, j, drop=drop]
343      else if(missing(j))      else if(missing(j))
344          x[i, , drop=drop]          if(nargs() == 4) x[i, , drop=drop] else x[i, drop=drop]
345      else      else
346          x[i,j, drop=drop]          x[i,j, drop=drop]
347      if(isS4(x) && isDiagonal(x)) as(x, "diagonalMatrix") else x      if(isS4(x) && isDiagonal(x)) as(x, "diagonalMatrix") else x
# Line 254  Line 351 
351                           j = "index", drop = "logical"), subDiag)                           j = "index", drop = "logical"), subDiag)
352  setMethod("[", signature(x = "diagonalMatrix", i = "index",  setMethod("[", signature(x = "diagonalMatrix", i = "index",
353                          j = "missing", drop = "logical"),                          j = "missing", drop = "logical"),
354            function(x, i, j, ..., drop) subDiag(x, i=i, drop=drop))            function(x, i, j, ..., drop) {
355                  na <- nargs()
356                  Matrix.msg("diag[i,m,l] : nargs()=", na, .M.level = 2)
357                  if(na == 4)
358                       subDiag(x, i=i, , drop=drop)
359                  else subDiag(x, i=i,   drop=drop)
360              })
361  setMethod("[", signature(x = "diagonalMatrix", i = "missing",  setMethod("[", signature(x = "diagonalMatrix", i = "missing",
362                           j = "index", drop = "logical"),                           j = "index", drop = "logical"),
363            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 367 
367  ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch  ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch
368  ## Only(?) current bug:  x[i] <- value  is wrong when  i is *vector*  ## Only(?) current bug:  x[i] <- value  is wrong when  i is *vector*
369  replDiag <- function(x, i, j, ..., value) {  replDiag <- function(x, i, j, ..., value) {
370      x <- as(x, "sparseMatrix")      x <- as(x, "TsparseMatrix")
371      if(missing(i))      if(missing(i))
372          x[, j] <- value          x[, j] <- value
373      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 412 
412                               x@x[ii] <- value                               x@x[ii] <- value
413                               x                               x
414                           } else { ## no longer diagonal, but remain sparse:                           } else { ## no longer diagonal, but remain sparse:
415                               x <- as(x, "sparseMatrix")                               x <- as(x, "TsparseMatrix")
416                               x[i] <- value                               x[i] <- value
417                               x                               x
418                           }                           }
# Line 325  Line 428 
428                                  j = "index", value = "replValue"),                                  j = "index", value = "replValue"),
429                   function(x,i,j, ..., value) replDiag(x, j=j, value=value))                   function(x,i,j, ..., value) replDiag(x, j=j, value=value))
430    
431    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",
432                                    value = "sparseMatrix"),
433                     function (x, i, j, ..., value)
434                     callGeneric(x=x, , j=j, value = as(value, "sparseVector")))
435    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "missing",
436                                    value = "sparseMatrix"),
437                     function (x, i, j, ..., value)
438                     callGeneric(x=x, i=i, , value = as(value, "sparseVector")))
439    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "index",
440                                    value = "sparseMatrix"),
441                     function (x, i, j, ..., value)
442                     callGeneric(x=x, i=i, j=j, value = as(value, "sparseVector")))
443    
444    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",
445                                    value = "sparseVector"),
446                     replDiag)
447    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "missing",
448                                    value = "sparseVector"),
449                     replDiag)
450    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "index",
451                                    value = "sparseVector"),
452                     replDiag)
453    
454    
455  setMethod("t", signature(x = "diagonalMatrix"),  setMethod("t", signature(x = "diagonalMatrix"),
456            function(x) { x@Dimnames <- x@Dimnames[2:1] ; x })            function(x) { x@Dimnames <- x@Dimnames[2:1] ; x })
# Line 472  Line 598 
598  ##        function(x, y = NULL) {  ##        function(x, y = NULL) {
599  ##           })  ##           })
600    
601  ## setMethod("tcrossprod", signature(x = "denseMatrix", y = "diagonalMatrix"),  Cspdiagprod <- function(x, y) {
602  ##        function(x, y = NULL) {      dx <- dim(x <- .Call(Csparse_diagU2N, x))
603  ##           })      dy <- dim(y)
604        if(dx[2] != dy[1]) stop("non-matching dimensions")
605        ind <- rep.int(seq_len(dx[2]), x@p[-1] - x@p[-dx[2]-1L])
606        if(y@diag == "N")
607            x@x <- x@x * y@x[ind]
608        x
609    }
610    
611    diagCspprod <- function(x, y) {
612        dx <- dim(x)
613        dy <- dim(y <- .Call(Csparse_diagU2N, y))
614        if(dx[2] != dy[1]) stop("non-matching dimensions")
615        if(x@diag == "N")
616            y@x <- y@x * x@x[y@i + 1L]
617        y
618    }
619    
620    setMethod("crossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
621              function(x, y = NULL) diagCspprod(x, y))
622    
623  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
624            function(x, y = NULL) crossprod(as(x, "sparseMatrix"), y))            function(x, y = NULL) diagCspprod(x, as(y, "CsparseMatrix")))
625    
626    ## Prefer calling diagCspprod to Cspdiagprod if going to transpose anyway
627    ##  x'y == (y'x)'
628    setMethod("crossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
629              function(x, y = NULL) t(diagCspprod(y, x)))
630    
631  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
632            function(x, y = NULL) crossprod(x, as(y, "sparseMatrix")))            function(x, y = NULL) t(diagCspprod(y, as(x, "Csparsematrix"))))
633    
634    setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
635              function(x, y = NULL) diagCspprod(x, t(y)))
636    
637  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
638            function(x, y = NULL) tcrossprod(as(x, "sparseMatrix"), y))            function(x, y = NULL) diagCspprod(x, t(as(y, "CsparseMatrix"))))
639    
640    setMethod("tcrossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
641              function(x, y = NULL) Cspdiagprod(x, y))
642    
643  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
644            function(x, y = NULL) tcrossprod(x, as(y, "sparseMatrix")))            function(x, y = NULL) Cspdiagprod(as(x, "CsparseMatrix"), y))
645    
646    setMethod("%*%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
647              function(x, y) diagCspprod(x, y))
648    
 ## FIXME?: In theory, this can be done *FASTER*, in some cases, via tapply1()  
649  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
650            function(x, y) as(x, "sparseMatrix") %*% y)            function(x, y) diagCspprod(as(x, "CsparseMatrix"), y))
651    
652  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
653            function(x, y) x %*% as(y, "sparseMatrix"))            function(x, y) Cspdiagprod(as(x, "CsparseMatrix"), y))
654  ## 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)  
655  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
656            function(x, y) x %*% as(y, "CsparseMatrix"))            function(x, y) Cspdiagprod(x, y))
657  ## NB: this is *not* needed for Tsparse & Rsparse  
658  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*
659  ##       do indeed work by going through sparse (and *not* ddense)!  ##       do indeed work by going through sparse (and *not* ddense)!
660    
   
   
661  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),
662            function(a, b, ...) {            function(a, b, ...) {
663                a@x <- 1/ a@x                a@x <- 1/ a@x
# Line 559  Line 710 
710          isNum <- (is.numeric(r) || is.numeric(r00))          isNum <- (is.numeric(r) || is.numeric(r00))
711          isLog <- (is.logical(r) || is.logical(r00))          isLog <- (is.logical(r) || is.logical(r00))
712    
713          if(getOption("verbose"))          Matrix.msg("exploding   <diag>  o  <diag>  into dense matrix")
             message("exploding  <diag>  o  <diag>  into dense matrix")  
714          d <- e1@Dim          d <- e1@Dim
715          n <- d[1]          n <- d[1]
716          stopifnot(length(r) == n)          stopifnot(length(r) == n)
# Line 619  Line 769 
769  ## ddi*: Arith: result numeric, potentially ddiMatrix  ## ddi*: Arith: result numeric, potentially ddiMatrix
770  setMethod("Arith", signature(e1 = "ddiMatrix", e2 = "numeric"),  setMethod("Arith", signature(e1 = "ddiMatrix", e2 = "numeric"),
771            function(e1,e2) {            function(e1,e2) {
772                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]; nsq <- n^2
773                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
774                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
775                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
# Line 637  Line 787 
787    
788  setMethod("Arith", signature(e1 = "numeric", e2 = "ddiMatrix"),  setMethod("Arith", signature(e1 = "numeric", e2 = "ddiMatrix"),
789            function(e1,e2) {            function(e1,e2) {
790                n <- e2@Dim[1]; nsq <- n*n                n <- e2@Dim[1]; nsq <- n^2
791                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
792                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
793                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
# Line 656  Line 806 
806  ## ldi* Arith --> result numeric, potentially ddiMatrix  ## ldi* Arith --> result numeric, potentially ddiMatrix
807  setMethod("Arith", signature(e1 = "ldiMatrix", e2 = "numeric"),  setMethod("Arith", signature(e1 = "ldiMatrix", e2 = "numeric"),
808            function(e1,e2) {            function(e1,e2) {
809                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]; nsq <- n^2
810                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
811                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
812                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
# Line 675  Line 825 
825    
826  setMethod("Arith", signature(e1 = "numeric", e2 = "ldiMatrix"),  setMethod("Arith", signature(e1 = "numeric", e2 = "ldiMatrix"),
827            function(e1,e2) {            function(e1,e2) {
828                n <- e2@Dim[1]; nsq <- n*n                n <- e2@Dim[1]; nsq <- n^2
829                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
830                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
831                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
# Line 695  Line 845 
845  ## ddi*: for "Ops" without Arith --> result logical, potentially ldi  ## ddi*: for "Ops" without Arith --> result logical, potentially ldi
846  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "numeric"),  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "numeric"),
847            function(e1,e2) {            function(e1,e2) {
848                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]; nsq <- n^2
849                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
850                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
851                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
# Line 714  Line 864 
864    
865  setMethod("Ops", signature(e1 = "numeric", e2 = "ddiMatrix"),  setMethod("Ops", signature(e1 = "numeric", e2 = "ddiMatrix"),
866            function(e1,e2) {            function(e1,e2) {
867                n <- e2@Dim[1]; nsq <- n*n                n <- e2@Dim[1]; nsq <- n^2
868                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
869                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
870                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
# Line 734  Line 884 
884  ## ldi*: for "Ops" without Arith --> result logical, potentially ldi  ## ldi*: for "Ops" without Arith --> result logical, potentially ldi
885  setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "numeric"),  setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "numeric"),
886            function(e1,e2) {            function(e1,e2) {
887                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]; nsq <- n^2
888                f0 <- callGeneric(FALSE, e2)                f0 <- callGeneric(FALSE, e2)
889                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
890                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
# Line 752  Line 902 
902    
903  setMethod("Ops", signature(e1 = "numeric", e2 = "ldiMatrix"),  setMethod("Ops", signature(e1 = "numeric", e2 = "ldiMatrix"),
904            function(e1,e2) {            function(e1,e2) {
905                n <- e2@Dim[1]; nsq <- n*n                n <- e2@Dim[1]; nsq <- n^2
906                f0 <- callGeneric(e1, FALSE)                f0 <- callGeneric(e1, FALSE)
907                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
908                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
# Line 783  Line 933 
933      setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),      setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),
934                function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "l")))                function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "l")))
935  }  }
936  ## This should *not* dispatch to <dense> methods (in ./Ops.R ), as  
937  ##  FALSE & <anything> |-> FALSE : hence result should be diagonal:  ## Direct subclasses of "denseMatrix": currently ddenseMatrix, ldense... :
938  for(cl in diCls) {  dense.subCl <- local({ dM.scl <- getClass("denseMatrix")@subclasses
939      setMethod("&", signature(e1 = cl, e2 = "ANY"),                         names(dM.scl)[sapply(dM.scl, slot, "distance") == 1] })
940                function(e1,e2) e1 & as(e2,"Matrix"))  for(DI in diCls) {
941      setMethod("&", signature(e1 = "ANY", e2 = cl),      for(c2 in c(dense.subCl, "Matrix")) {
942                function(e1,e2) as(e1,"Matrix") & e2)          for(Fun in c("*", "^", "&")) {
943      for(c2 in c("denseMatrix", "Matrix")) {              setMethod(Fun, signature(e1 = DI, e2 = c2),
944          setMethod("&", signature(e1 = cl, e2 = c2),                        function(e1,e2) callGeneric(e1, Diagonal(x = diag(e2))))
945                    function(e1,e2) e1 & Diagonal(x = diag(e2)))              setMethod(Fun, signature(e1 = c2, e2 = DI),
946          setMethod("&", signature(e1 = c2, e2 = cl),                        function(e1,e2) callGeneric(Diagonal(x = diag(e1)), e2))
947                    function(e1,e2) Diagonal(x = diag(e1)) & e2)          }
948            ## NB: This arguably implicitly uses  0/0 :== 0  to keep diagonality
949            for(Fun in c("%%", "%/%", "/")) {
950                setMethod(Fun, signature(e1 = DI, e2 = c2),
951                          function(e1,e2) callGeneric(e1, Diagonal(x = diag(e2))))
952            }
953      }      }
954  }  }
955    
956    
957  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"
958  ### ----------  any, all: separately here  ### ----------   the last 4: separately here
959  for(cl in diCls) {  for(cl in diCls) {
960  setMethod("any", cl,  setMethod("any", cl,
961            function (x, ..., na.rm) {            function (x, ..., na.rm) {
962                if(any(x@Dim == 0)) FALSE                if(any(x@Dim == 0)) FALSE
963                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)
964            })            })
965  setMethod("all",  cl, function (x, ..., na.rm) any(x@Dim == 0))  setMethod("all",  cl, function (x, ..., na.rm) {
966  setMethod("prod", cl, function (x, ..., na.rm) as.numeric(any(x@Dim == 0)))      n <- x@Dim[1]
967        if(n >= 2) FALSE
968        else if(n == 0 || x@diag == "U") TRUE
969        else all(x@x, ..., na.rm = na.rm)
970    })
971    setMethod("prod", cl, function (x, ..., na.rm) {
972        n <- x@Dim[1]
973        if(n >= 2) 0
974        else if(n == 0 || x@diag == "U") 1
975        else ## n == 1, diag = "N" :
976            prod(x@x, ..., na.rm = na.rm)
977    })
978    
979  setMethod("sum", cl,  setMethod("sum", cl,
980            function(x, ..., na.rm) {            function(x, ..., na.rm) {
# Line 847  Line 1013 
1013      invisible(x)      invisible(x)
1014  }  }
1015    
1016    ## somewhat consistent with "print" for sparseMatrix :
1017    setMethod("print", signature(x = "diagonalMatrix"), prDiag)
1018    
1019  setMethod("show", signature(object = "diagonalMatrix"),  setMethod("show", signature(object = "diagonalMatrix"),
1020            function(object) {            function(object) {
1021                d <- dim(object)                d <- dim(object)
1022                cl <- class(object)                cl <- class(object)
1023                cat(sprintf('%d x %d diagonal matrix of class "%s"\n',                cat(sprintf('%d x %d diagonal matrix of class "%s"',
1024                            d[1], d[2], cl))                            d[1], d[2], cl))
1025                  if(d[1] < 50) {
1026                      cat("\n")
1027                prDiag(object)                prDiag(object)
1028                  } else {
1029                      cat(", with diagonal entries\n")
1030                      show(diag(object))
1031                      invisible(object)
1032                  }
1033            })            })

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

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