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

pkg/R/diagMatrix.R revision 2197, Mon Jun 2 14:34:42 2008 UTC pkg/Matrix/R/diagMatrix.R revision 2741, Fri Dec 9 10:38:51 2011 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,m), uplo = "U",
36  .symDiagonal <- function(n, x = rep.int(1,n), uplo = "U") {                              shape = if(missing(cols)) "t" else "g",
37                                kind, cols = if(n) 0:(n - 1L) else integer(0))
38    {
39      stopifnot(n == (n. <- as.integer(n)), (n <- n.) >= 0)      stopifnot(n == (n. <- as.integer(n)), (n <- n.) >= 0)
40      if((lx <- length(x)) == 1) x <- rep.int(x, n)      if(!missing(cols))
41      else if(lx != n) stop("length(x) must be 1 or n")          stopifnot(0 <= (cols <- as.integer(cols)), cols < n)
42      cls <-      m <- length(cols)
43          if(is.double(x)) "dsCMatrix"      if(missing(kind))
44          else if(is.logical(x)) "lsCMatrix"          kind <-
45                if(is.double(x)) "d"
46                else if(is.logical(x)) "l"
47          else { ## for now          else { ## for now
48              storage.mode(x) <- "double"              storage.mode(x) <- "double"
49              "dsCMatrix"                  "d"
50          }          }
51      new(cls, Dim = c(n,n), x = x, uplo = uplo,      else stopifnot(any(kind == c("d","l","n")))
52          i = if(n) 0:(n - 1L) else integer(0), p = 0:n)      if(kind != "n") {
53            if((lx <- length(x)) == 1) x <- rep.int(x, m)
54            else if(lx != m) stop("length(x) must be either 1 or #{cols}")
55        }
56        stopifnot(is.character(shape), nchar(shape) == 1,
57                  any(shape == c("t","s","g"))) # triangular / symmetric / general
58        if(kind == "n") {
59            if(shape == "g")
60                new("ngCMatrix", Dim = c(n,m), i = cols, p = 0:m)
61            else new(paste0("n", shape, "CMatrix"), Dim = c(n,m), uplo = uplo,
62                     i = cols, p = 0:m)
63        }
64        ## kind != "n" -- have x slot :
65        else if(shape == "g")
66            new(paste0(kind, "gCMatrix"), Dim = c(n,m),
67                x = x, i = cols, p = 0:m)
68        else new(paste0(kind, shape, "CMatrix"), Dim = c(n,m), uplo = uplo,
69                 x = x, i = cols, p = 0:m)
70  }  }
71    
72    ## Pkg 'spdep' had (relatively slow) versions of this as_dsCMatrix_I()
73    .symDiagonal <- function(n, x = rep.int(1,n), uplo = "U")
74        .sparseDiagonal(n, x, uplo, shape = "s")
75    
76    # instead of   diagU2N(as(Diagonal(n), "CsparseMatrix")), diag = "N" in any case:
77    .trDiagonal <- function(n, x = rep.int(1,n), uplo = "U")
78        .sparseDiagonal(n, x, uplo, shape = "t")
79    
80    
81  ### 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.
82  ### 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
83  ### 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
84  ### who posted his bdiag() function written in December 1995.  ### who posted his bdiag() function written in December 1995.
85    if(FALSE)##--- no longer used:
86  bdiag <- function(...) {  .bdiag <- function(lst) {
87      if(nargs() == 0) return(new("dgCMatrix"))      ### block-diagonal matrix [a dgTMatrix] from list of matrices
88      ## else :      stopifnot(is.list(lst), length(lst) >= 1)
89      mlist <- if (nargs() == 1) as.list(...) else list(...)      dims <- sapply(lst, dim, USE.NAMES=FALSE)
     dims <- sapply(mlist, dim)  
90      ## make sure we had all matrices:      ## make sure we had all matrices:
91      if(!(is.matrix(dims) && nrow(dims) == 2))      if(!(is.matrix(dims) && nrow(dims) == 2))
92          stop("some arguments are not matrices")          stop("some arguments are not matrices")
93      csdim <- rbind(rep.int(0L, 2),      csdim <- rbind(rep.int(0L, 2),
94                     apply(sapply(mlist, dim), 1, cumsum))                     apply(dims, 1, cumsum))
95      ret <- new("dgTMatrix", Dim = as.integer(csdim[nrow(csdim),]))      r <- new("dgTMatrix")
96        r@Dim <- as.integer(csdim[nrow(csdim),])
97      add1 <- matrix(1:0, 2,2)      add1 <- matrix(1:0, 2,2)
98      for(i in seq_along(mlist)) {      for(i in seq_along(lst)) {
99          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])
100          if(is.null(dim(indx))) ## non-square matrix          if(is.null(dim(indx))) ## non-square matrix
101              ret[indx[[1]],indx[[2]]] <- mlist[[i]]              r[indx[[1]],indx[[2]]] <- lst[[i]]
102          else ## square matrix          else ## square matrix
103              ret[indx[,1],indx[,2]] <- mlist[[i]]              r[indx[,1], indx[,2]] <- lst[[i]]
104        }
105        r
106    }
107    ## expand(<mer>) needed something like bdiag() for lower-triangular
108    ## (Tsparse) Matrices; hence Doug Bates provided a much more efficient
109    ##  implementation for those; now extended and generalized:
110    .bdiag <- function(lst) {
111        ## block-diagonal matrix [a dgTMatrix] from list of matrices
112        stopifnot(is.list(lst), (nl <- length(lst)) >= 1)
113    
114        Tlst <- lapply(lapply(lst, Matrix:::as_Csp2), # includes "diagU2N"
115                       as, "TsparseMatrix")
116        if(nl == 1) return(Tlst[[1]])
117        ## else
118        i_off <- c(0L, cumsum(sapply(Tlst, nrow)))
119        j_off <- c(0L, cumsum(sapply(Tlst, ncol)))
120    
121        clss <- sapply(Tlst, class)
122        typ <- substr(clss, 2, 2)
123        knd <- substr(clss, 1, 1)
124        sym <- typ == "s" # symmetric ones
125        tri <- typ == "t" # triangular ones
126        use.n <- any(is.n <- knd == "n")
127        if(use.n && !(use.n <- all(is.n))) {
128            Tlst[is.n] <- lapply(Tlst[is.n], as, "lMatrix")
129            knd [is.n] <- "l"
130        }
131        use.l <- !use.n && all(knd == "l")
132        if(all(sym)) { ## result should be *symmetric*
133            uplos <- sapply(Tlst, slot, "uplo") ## either "U" or "L"
134            tLU <- table(uplos)# of length 1 or 2 ..
135            if(length(tLU) == 1) { ## all "U" or all "L"
136                useU <- uplos[1] == "U"
137            } else { ## length(tLU) == 2, counting "L" and "U"
138                useU <- diff(tLU) >= 0
139                if(useU && (hasL <- tLU[1] > 0))
140                    Tlst[hasL] <- lapply(Tlst[hasL], t)
141                else if(!useU && (hasU <- tLU[2] > 0))
142                    Tlst[hasU] <- lapply(Tlst[hasU], t)
143            }
144            if(use.n) { ## return nsparseMatrix :
145                r <- new("nsTMatrix")
146            } else {
147                r <- new(paste0(if(use.l) "l" else "d", "sTMatrix"))
148                r@x <- unlist(lapply(Tlst, slot, "x"))
149            }
150            r@uplo <- if(useU) "U" else "L"
151        }
152        else if(all(tri) && { ULs <- sapply(Tlst, slot, "uplo")##  "U" or "L"
153                              all(ULs[1L] == ULs[-1L]) } ## all upper or all lower
154           ){ ## *triangular* result
155    
156            if(use.n) { ## return nsparseMatrix :
157                r <- new("ntTMatrix")
158            } else {
159                r <- new(paste0(if(use.l) "l" else "d", "tTMatrix"))
160                r@x <- unlist(lapply(Tlst, slot, "x"))
161            }
162            r@uplo <- ULs[1L]
163        }
164        else {
165            if(any(sym))
166                Tlst[sym] <- lapply(Tlst[sym], as, "generalMatrix")
167            if(use.n) { ## return nsparseMatrix :
168                r <- new("ngTMatrix")
169            } else {
170                r <- new(paste0(if(use.l) "l" else "d", "gTMatrix"))
171                r@x <- unlist(lapply(Tlst, slot, "x"))
172            }
173      }      }
174      ## slightly debatable if we really should return Csparse.. :      r@Dim <- c(i_off[nl+1], j_off[nl + 1])
175      as(ret, "CsparseMatrix")      r@i <- unlist(lapply(1:nl, function(k) Tlst[[k]]@i + i_off[k]))
176        r@j <- unlist(lapply(1:nl, function(k) Tlst[[k]]@j + j_off[k]))
177        r
178    }
179    
180    bdiag <- function(...) {
181        if((nA <- nargs()) == 0) return(new("dgCMatrix"))
182        if(nA == 1 && !is.list(...))
183            return(as(..., "CsparseMatrix"))
184        alis <- if(nA == 1 && is.list(..1)) ..1 else list(...)
185        if(length(alis) == 1)
186            return(as(alis[[1]], "CsparseMatrix"))
187    
188        ## else : two or more arguments
189        as(.bdiag(alis), "CsparseMatrix")
190  }  }
191    
192    
# Line 124  Line 237 
237  ## ddi*:  ## ddi*:
238  diag2tT <- function(from) .diag2tT(from, "U", "d")  diag2tT <- function(from) .diag2tT(from, "U", "d")
239  setAs("ddiMatrix", "triangularMatrix", diag2tT)  setAs("ddiMatrix", "triangularMatrix", diag2tT)
240  setAs("ddiMatrix", "sparseMatrix", diag2tT)  ##_no_longer_ setAs("ddiMatrix", "sparseMatrix", diag2tT)
241  ## needed too (otherwise <dense> -> Tsparse is taken):  ## needed too (otherwise <dense> -> Tsparse is taken):
242  setAs("ddiMatrix", "TsparseMatrix", diag2tT)  setAs("ddiMatrix", "TsparseMatrix", diag2tT)
243    setAs("ddiMatrix", "dsparseMatrix", diag2tT)
244  setAs("ddiMatrix", "CsparseMatrix",  setAs("ddiMatrix", "CsparseMatrix",
245        function(from) as(.diag2tT(from, "U", "d"), "CsparseMatrix"))        function(from) as(.diag2tT(from, "U", "d"), "CsparseMatrix"))
246  setAs("ddiMatrix", "symmetricMatrix",  setAs("ddiMatrix", "symmetricMatrix",
# Line 135  Line 249 
249  ## ldi*:  ## ldi*:
250  diag2tT <- function(from) .diag2tT(from, "U", "l")  diag2tT <- function(from) .diag2tT(from, "U", "l")
251  setAs("ldiMatrix", "triangularMatrix", diag2tT)  setAs("ldiMatrix", "triangularMatrix", diag2tT)
252  setAs("ldiMatrix", "sparseMatrix", diag2tT)  ##_no_longer_ setAs("ldiMatrix", "sparseMatrix", diag2tT)
253  ## needed too (otherwise <dense> -> Tsparse is taken):  ## needed too (otherwise <dense> -> Tsparse is taken):
254  setAs("ldiMatrix", "TsparseMatrix", diag2tT)  setAs("ldiMatrix", "TsparseMatrix", diag2tT)
255    setAs("ldiMatrix", "lsparseMatrix", diag2tT)
256  setAs("ldiMatrix", "CsparseMatrix",  setAs("ldiMatrix", "CsparseMatrix",
257        function(from) as(.diag2tT(from, "U", "l"), "CsparseMatrix"))        function(from) as(.diag2tT(from, "U", "l"), "CsparseMatrix"))
258  setAs("ldiMatrix", "symmetricMatrix",  setAs("ldiMatrix", "symmetricMatrix",
# Line 152  Line 267 
267                Dim = from@Dim, Dimnames = from@Dimnames)                Dim = from@Dim, Dimnames = from@Dimnames)
268        })        })
269    
270    setAs("diagonalMatrix", "nsparseMatrix", function(from) as(from, "nMatrix"))
271    
272  ## Cheap fast substitute for diag() which *does* preserve the mode of x :  ## Cheap fast substitute for diag() which *does* preserve the mode of x :
273  mkDiag <- function(x, n) {  mkDiag <- function(x, n) {
# Line 181  Line 297 
297  setAs("diagonalMatrix", "generalMatrix", # prefer sparse:  setAs("diagonalMatrix", "generalMatrix", # prefer sparse:
298        function(from) as(as(from, "CsparseMatrix"), "generalMatrix"))        function(from) as(as(from, "CsparseMatrix"), "generalMatrix"))
299    
300    setAs("diagonalMatrix", "denseMatrix",
301          function(from) as(as(from, "CsparseMatrix"), "denseMatrix"))
302    
303  .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
304    
305  .diag.2N <- function(m) {  .diag.2N <- function(m) {
# Line 201  Line 320 
320            d <- dim(from)            d <- dim(from)
321            if(d[1] != (n <- d[2])) stop("non-square matrix")            if(d[1] != (n <- d[2])) stop("non-square matrix")
322            if(any(from[row(from) != col(from)] != 0))            if(any(from[row(from) != col(from)] != 0))
323                stop("has non-zero off-diagonal entries")                stop("matrix with non-zero off-diagonals cannot be coerced to diagonalMatrix")
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) ## uni := {is it unit-diagonal ?}
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            } ## 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",
# Line 225  Line 344 
344            x <- diag(from)            x <- diag(from)
345            if(is.logical(x)) {            if(is.logical(x)) {
346                cl <- "ldiMatrix"                cl <- "ldiMatrix"
347                uni <- all(x)                uni <- allTrue(x)
348            } else {            } else {
349                cl <- "ddiMatrix"                cl <- "ddiMatrix"
350                uni <- all(x == 1)                uni <- allTrue(x == 1)
351                storage.mode(x) <- "double"                storage.mode(x) <- "double"
352            }            } ## TODO: complex
353            new(cl, Dim = c(n,n), diag = if(uni) "U" else "N",            new(cl, Dim = c(n,n), diag = if(uni) "U" else "N",
354                x = if(uni) x[FALSE] else x)                x = if(uni) x[FALSE] else x)
355        })        })
# Line 240  Line 359 
359            function(x = 1, nrow, ncol) .diag.x(x))            function(x = 1, nrow, ncol) .diag.x(x))
360    
361  subDiag <- function(x, i, j, ..., drop) {  subDiag <- function(x, i, j, ..., drop) {
362      x <- as(x, "sparseMatrix")      x <- as(x, "CsparseMatrix") ## << was "TsparseMatrix" (Csparse is faster now)
363      x <- if(missing(i))      x <- if(missing(i))
364          x[, j, drop=drop]          x[, j, drop=drop]
365      else if(missing(j))      else if(missing(j))
366          x[i, , drop=drop]          if(nargs() == 4) x[i, , drop=drop] else x[i, drop=drop]
367      else      else
368          x[i,j, drop=drop]          x[i,j, drop=drop]
369      if(isS4(x) && isDiagonal(x)) as(x, "diagonalMatrix") else x      if(isS4(x) && isDiagonal(x)) as(x, "diagonalMatrix") else x
# Line 254  Line 373 
373                           j = "index", drop = "logical"), subDiag)                           j = "index", drop = "logical"), subDiag)
374  setMethod("[", signature(x = "diagonalMatrix", i = "index",  setMethod("[", signature(x = "diagonalMatrix", i = "index",
375                          j = "missing", drop = "logical"),                          j = "missing", drop = "logical"),
376            function(x, i, j, ..., drop) subDiag(x, i=i, drop=drop))            function(x, i, j, ..., drop) {
377                  na <- nargs()
378                  Matrix.msg("diag[i,m,l] : nargs()=", na, .M.level = 2)
379                  if(na == 4)
380                       subDiag(x, i=i, , drop=drop)
381                  else subDiag(x, i=i,   drop=drop)
382              })
383  setMethod("[", signature(x = "diagonalMatrix", i = "missing",  setMethod("[", signature(x = "diagonalMatrix", i = "missing",
384                           j = "index", drop = "logical"),                           j = "index", drop = "logical"),
385            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 389 
389  ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch  ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch
390  ## Only(?) current bug:  x[i] <- value  is wrong when  i is *vector*  ## Only(?) current bug:  x[i] <- value  is wrong when  i is *vector*
391  replDiag <- function(x, i, j, ..., value) {  replDiag <- function(x, i, j, ..., value) {
392      x <- as(x, "sparseMatrix")      x <- as(x, "TsparseMatrix")
393      if(missing(i))      if(missing(i))
394          x[, j] <- value          x[, j] <- value
395      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 434 
434                               x@x[ii] <- value                               x@x[ii] <- value
435                               x                               x
436                           } else { ## no longer diagonal, but remain sparse:                           } else { ## no longer diagonal, but remain sparse:
437                               x <- as(x, "sparseMatrix")                               x <- as(x, "TsparseMatrix")
438                               x[i] <- value                               x[i] <- value
439                               x                               x
440                           }                           }
# Line 325  Line 450 
450                                  j = "index", value = "replValue"),                                  j = "index", value = "replValue"),
451                   function(x,i,j, ..., value) replDiag(x, j=j, value=value))                   function(x,i,j, ..., value) replDiag(x, j=j, value=value))
452    
453    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",
454                                    value = "sparseMatrix"),
455                     function (x, i, j, ..., value)
456                     callGeneric(x=x, , j=j, value = as(value, "sparseVector")))
457    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "missing",
458                                    value = "sparseMatrix"),
459                     function (x, i, j, ..., value)
460                     callGeneric(x=x, i=i, , value = as(value, "sparseVector")))
461    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "index",
462                                    value = "sparseMatrix"),
463                     function (x, i, j, ..., value)
464                     callGeneric(x=x, i=i, j=j, value = as(value, "sparseVector")))
465    
466    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",
467                                    value = "sparseVector"),
468                     replDiag)
469    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "missing",
470                                    value = "sparseVector"),
471                     replDiag)
472    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "index",
473                                    value = "sparseVector"),
474                     replDiag)
475    
476    
477  setMethod("t", signature(x = "diagonalMatrix"),  setMethod("t", signature(x = "diagonalMatrix"),
478            function(x) { x@Dimnames <- x@Dimnames[2:1] ; x })            function(x) { x@Dimnames <- x@Dimnames[2:1] ; x })
# Line 371  Line 519 
519  ##       ---------------------  ##       ---------------------
520  ## Note that "ldi" logical are treated as numeric  ## Note that "ldi" logical are treated as numeric
521  diagdiagprod <- function(x, y) {  diagdiagprod <- function(x, y) {
522      n <- dimCheck(x,y)[1]      dimCheck(x,y)
523      if(x@diag != "U") {      if(x@diag != "U") {
524          if(y@diag != "U") {          if(y@diag != "U") {
525              nx <- x@x * y@x              nx <- x@x * y@x
# Line 403  Line 551 
551      dx <- dim(x)      dx <- dim(x)
552      dy <- dim(y)      dy <- dim(y)
553      if(dx[2] != dy[1]) stop("non-matching dimensions")      if(dx[2] != dy[1]) stop("non-matching dimensions")
     n <- dx[1]  
554      as(if(x@diag == "U") y else x@x * y, "Matrix")      as(if(x@diag == "U") y else x@x * y, "Matrix")
555  }  }
556  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),
# Line 472  Line 619 
619  ##        function(x, y = NULL) {  ##        function(x, y = NULL) {
620  ##           })  ##           })
621    
622  ## setMethod("tcrossprod", signature(x = "denseMatrix", y = "diagonalMatrix"),  Cspdiagprod <- function(x, y) {
623  ##        function(x, y = NULL) {      dx <- dim(x <- .Call(Csparse_diagU2N, x))
624  ##           })      dy <- dim(y)
625        if(dx[2] != dy[1]) stop("non-matching dimensions")
626        if(y@diag == "N") {
627            if(!all(y@x[1L] == y@x[-1L]) && is(x, "symmetricMatrix"))
628                x <- as(x, "generalMatrix")
629            ind <- rep.int(seq_len(dx[2]), x@p[-1] - x@p[-dx[2]-1L])
630            x@x <- x@x * y@x[ind]
631        }
632        if(is(x, "compMatrix") && length(xf <- x@factors)) {
633            ## instead of dropping all factors, be smart about some
634            ## TODO ......
635            x@factors <- list()
636        }
637        x
638    }
639    
640    diagCspprod <- function(x, y) {
641        dx <- dim(x)
642        dy <- dim(y <- .Call(Csparse_diagU2N, y))
643        if(dx[2] != dy[1]) stop("non-matching dimensions")
644        if(x@diag == "N") {
645            if(!all(x@x[1L] == x@x[-1L]) && is(y, "symmetricMatrix"))
646                y <- as(y, "generalMatrix")
647            y@x <- y@x * x@x[y@i + 1L]
648        }
649        if(is(y, "compMatrix") && length(yf <- y@factors)) {
650            ## instead of dropping all factors, be smart about some
651            ## TODO
652            keep <- character()
653            if(iLU <- names(yf) == "LU") {
654                ## TODO keep <- "LU"
655            }
656            y@factors <- yf[keep]
657        }
658        y
659    }
660    
661    setMethod("crossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
662              function(x, y = NULL) diagCspprod(x, y))
663    
664  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
665            function(x, y = NULL) crossprod(as(x, "sparseMatrix"), y))            function(x, y = NULL) diagCspprod(x, as(y, "CsparseMatrix")))
666    
667    ## Prefer calling diagCspprod to Cspdiagprod if going to transpose anyway
668    ##  x'y == (y'x)'
669    setMethod("crossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
670              function(x, y = NULL) t(diagCspprod(y, x)))
671    
672  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
673            function(x, y = NULL) crossprod(x, as(y, "sparseMatrix")))            function(x, y = NULL) t(diagCspprod(y, as(x, "Csparsematrix"))))
674    
675    setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
676              function(x, y = NULL) diagCspprod(x, t(y)))
677    
678  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
679            function(x, y = NULL) tcrossprod(as(x, "sparseMatrix"), y))            function(x, y = NULL) diagCspprod(x, t(as(y, "CsparseMatrix"))))
680    
681    setMethod("tcrossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
682              function(x, y = NULL) Cspdiagprod(x, y))
683    
684  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
685            function(x, y = NULL) tcrossprod(x, as(y, "sparseMatrix")))            function(x, y = NULL) Cspdiagprod(as(x, "CsparseMatrix"), y))
686    
687    setMethod("%*%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
688              function(x, y) diagCspprod(x, y))
689    
 ## FIXME?: In theory, this can be done *FASTER*, in some cases, via tapply1()  
690  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
691            function(x, y) as(x, "sparseMatrix") %*% y)            function(x, y) diagCspprod(as(x, "CsparseMatrix"), y))
692    
693  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
694            function(x, y) x %*% as(y, "sparseMatrix"))            function(x, y) Cspdiagprod(as(x, "CsparseMatrix"), y))
695  ## 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)  
696  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
697            function(x, y) x %*% as(y, "CsparseMatrix"))            function(x, y) Cspdiagprod(x, y))
698  ## NB: this is *not* needed for Tsparse & Rsparse  
699  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*
700  ##       do indeed work by going through sparse (and *not* ddense)!  ##       do indeed work by going through sparse (and *not* ddense)!
701    
   
   
702  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),
703            function(a, b, ...) {            function(a, b, ...) {
704                a@x <- 1/ a@x                a@x <- 1/ a@x
# Line 515  Line 707 
707            })            })
708    
709  solveDiag <- function(a, b, ...) {  solveDiag <- function(a, b, ...) {
710      if((n <- a@Dim[1]) != nrow(b))      if(a@Dim[1] != nrow(b))
711          stop("incompatible matrix dimensions")          stop("incompatible matrix dimensions")
712      ## trivially invert a 'in place' and multiply:      ## trivially invert a 'in place' and multiply:
713      a@x <- 1/ a@x      a@x <- 1/ a@x
# Line 556  Line 748 
748          .diag.2N(e1)          .diag.2N(e1)
749      }      }
750      else { ## result not diagonal, but at least symmetric:      else { ## result not diagonal, but at least symmetric:
751            ## e.g., m == m
752          isNum <- (is.numeric(r) || is.numeric(r00))          isNum <- (is.numeric(r) || is.numeric(r00))
753          isLog <- (is.logical(r) || is.logical(r00))          isLog <- (is.logical(r) || is.logical(r00))
754            Matrix.msg("exploding <diag> o <diag> into dense matrix", .M.level = 2)
         if(getOption("verbose"))  
             message("exploding  <diag>  o  <diag>  into dense matrix")  
755          d <- e1@Dim          d <- e1@Dim
756          n <- d[1]          n <- d[1]
757          stopifnot(length(r) == n)          stopifnot(length(r) == n)
# Line 588  Line 779 
779              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)
780  }  }
781    
782  ## FIXME:    diagonal  o  triangular  |-->  triangular  ## diagonal  o  triangular  |-->  triangular
783  ## -----     diagonal  o  symmetric   |-->  symmetric  ## diagonal  o  symmetric   |-->  symmetric
784  ##    {also when other is sparse: do these "here" --  ##    {also when other is sparse: do these "here" --
785  ##     before conversion to sparse, since that loses "diagonality"}  ##     before conversion to sparse, since that loses "diagonality"}
786    
# Line 619  Line 810 
810  ## ddi*: Arith: result numeric, potentially ddiMatrix  ## ddi*: Arith: result numeric, potentially ddiMatrix
811  setMethod("Arith", signature(e1 = "ddiMatrix", e2 = "numeric"),  setMethod("Arith", signature(e1 = "ddiMatrix", e2 = "numeric"),
812            function(e1,e2) {            function(e1,e2) {
813                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]; nsq <- n^2
814                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
815                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
816                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
# Line 637  Line 828 
828    
829  setMethod("Arith", signature(e1 = "numeric", e2 = "ddiMatrix"),  setMethod("Arith", signature(e1 = "numeric", e2 = "ddiMatrix"),
830            function(e1,e2) {            function(e1,e2) {
831                n <- e2@Dim[1]; nsq <- n*n                n <- e2@Dim[1]; nsq <- n^2
832                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
833                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
834                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
# Line 656  Line 847 
847  ## ldi* Arith --> result numeric, potentially ddiMatrix  ## ldi* Arith --> result numeric, potentially ddiMatrix
848  setMethod("Arith", signature(e1 = "ldiMatrix", e2 = "numeric"),  setMethod("Arith", signature(e1 = "ldiMatrix", 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 675  Line 866 
866    
867  setMethod("Arith", signature(e1 = "numeric", e2 = "ldiMatrix"),  setMethod("Arith", signature(e1 = "numeric", e2 = "ldiMatrix"),
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 695  Line 886 
886  ## ddi*: for "Ops" without Arith --> result logical, potentially ldi  ## ddi*: for "Ops" without Arith --> result logical, potentially ldi
887  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "numeric"),  setMethod("Ops", signature(e1 = "ddiMatrix", 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(0, e2)                f0 <- callGeneric(0, 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 714  Line 905 
905    
906  setMethod("Ops", signature(e1 = "numeric", e2 = "ddiMatrix"),  setMethod("Ops", signature(e1 = "numeric", e2 = "ddiMatrix"),
907            function(e1,e2) {            function(e1,e2) {
908                n <- e2@Dim[1]; nsq <- n*n                n <- e2@Dim[1]; nsq <- n^2
909                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
910                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
911                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
# Line 734  Line 925 
925  ## ldi*: for "Ops" without Arith --> result logical, potentially ldi  ## ldi*: for "Ops" without Arith --> result logical, potentially ldi
926  setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "numeric"),  setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "numeric"),
927            function(e1,e2) {            function(e1,e2) {
928                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]; nsq <- n^2
929                f0 <- callGeneric(FALSE, e2)                f0 <- callGeneric(FALSE, e2)
930                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
931                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
# Line 752  Line 943 
943    
944  setMethod("Ops", signature(e1 = "numeric", e2 = "ldiMatrix"),  setMethod("Ops", signature(e1 = "numeric", e2 = "ldiMatrix"),
945            function(e1,e2) {            function(e1,e2) {
946                n <- e2@Dim[1]; nsq <- n*n                n <- e2@Dim[1]; nsq <- n^2
947                f0 <- callGeneric(e1, FALSE)                f0 <- callGeneric(e1, FALSE)
948                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
949                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
# Line 774  Line 965 
965  for(other in c("ANY", "Matrix", "dMatrix")) {  for(other in c("ANY", "Matrix", "dMatrix")) {
966      ## ddi*:      ## ddi*:
967      setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other),      setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other),
968                function(e1,e2) callGeneric(diag2tT.u(e1,e2, "d"), e2))                function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2))
969      setMethod("Ops", signature(e1 = other, e2 = "ddiMatrix"),      setMethod("Ops", signature(e1 = other, e2 = "ddiMatrix"),
970                function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "d")))                function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "d")))
971      ## ldi*:      ## ldi*:
972      setMethod("Ops", signature(e1 = "ldiMatrix", e2 = other),      setMethod("Ops", signature(e1 = "ldiMatrix", e2 = other),
973                function(e1,e2) callGeneric(diag2tT.u(e1,e2, "l"), e2))                function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2))
974      setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),      setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),
975                function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "l")))                function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "l")))
976    }
977    
978    ## Direct subclasses of "denseMatrix": currently ddenseMatrix, ldense... :
979    if(FALSE)## too general, would contain  denseModelMatrix:
980    dense.subCl <- local({ dM.scl <- getClass("denseMatrix")@subclasses
981                           names(dM.scl)[sapply(dM.scl, slot, "distance") == 1] })
982    dense.subCl <- paste(c("d","l","n"), "denseMatrix", sep="")
983    for(DI in diCls) {
984        for(c2 in c(dense.subCl, "Matrix")) {
985            for(Fun in c("*", "^", "&")) {
986                setMethod(Fun, signature(e1 = DI, e2 = c2),
987                          function(e1,e2) callGeneric(e1, Diagonal(x = diag(e2))))
988                setMethod(Fun, signature(e1 = c2, e2 = DI),
989                          function(e1,e2) callGeneric(Diagonal(x = diag(e1)), e2))
990            }
991            ## NB: This arguably implicitly uses  0/0 :== 0  to keep diagonality
992            for(Fun in c("%%", "%/%", "/")) {
993                setMethod(Fun, signature(e1 = DI, e2 = c2),
994                          function(e1,e2) callGeneric(e1, Diagonal(x = diag(e2))))
995  }  }
 ## 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)  
996      }      }
997  }  }
998    
999    
1000  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"
1001  ### ----------  any, all: separately here  ### ----------   the last 4: separately here
1002  for(cl in diCls) {  for(cl in diCls) {
1003  setMethod("any", cl,  setMethod("any", cl,
1004            function (x, ..., na.rm) {            function (x, ..., na.rm) {
1005                if(any(x@Dim == 0)) FALSE                if(any(x@Dim == 0)) FALSE
1006                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)
1007            })            })
1008  setMethod("all",  cl, function (x, ..., na.rm) any(x@Dim == 0))  setMethod("all",  cl, function (x, ..., na.rm) {
1009  setMethod("prod", cl, function (x, ..., na.rm) as.numeric(any(x@Dim == 0)))      n <- x@Dim[1]
1010        if(n >= 2) FALSE
1011        else if(n == 0 || x@diag == "U") TRUE
1012        else all(x@x, ..., na.rm = na.rm)
1013    })
1014    setMethod("prod", cl, function (x, ..., na.rm) {
1015        n <- x@Dim[1]
1016        if(n >= 2) 0
1017        else if(n == 0 || x@diag == "U") 1
1018        else ## n == 1, diag = "N" :
1019            prod(x@x, ..., na.rm = na.rm)
1020    })
1021    
1022  setMethod("sum", cl,  setMethod("sum", cl,
1023            function(x, ..., na.rm) {            function(x, ..., na.rm) {
# Line 847  Line 1056 
1056      invisible(x)      invisible(x)
1057  }  }
1058    
1059    ## somewhat consistent with "print" for sparseMatrix :
1060    setMethod("print", signature(x = "diagonalMatrix"), prDiag)
1061    
1062  setMethod("show", signature(object = "diagonalMatrix"),  setMethod("show", signature(object = "diagonalMatrix"),
1063            function(object) {            function(object) {
1064                d <- dim(object)                d <- dim(object)
1065                cl <- class(object)                cl <- class(object)
1066                cat(sprintf('%d x %d diagonal matrix of class "%s"\n',                cat(sprintf('%d x %d diagonal matrix of class "%s"',
1067                            d[1], d[2], cl))                            d[1], d[2], cl))
1068                  if(d[1] < 50) {
1069                      cat("\n")
1070                prDiag(object)                prDiag(object)
1071                  } else {
1072                      cat(", with diagonal entries\n")
1073                      show(diag(object))
1074                      invisible(object)
1075                  }
1076            })            })

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

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