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 2633, Sun Dec 12 21:23:17 2010 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      ## slightly debatable if we really should return Csparse.. :      r
106      as(ret, "CsparseMatrix")  }
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        knds <- substr(clss, 2, 2)
123        sym  <- knds == "s" # symmetric ones
124        tri  <- knds == "t" # triangular ones
125        use.n <- any(is.n <- substr(clss,1,1) == "n")
126        if(use.n && !(use.n <- all(is.n)))
127            Tlst[is.n] <- lapply(Tlst[is.n], as, "lMatrix")
128        if(all(sym)) { ## result should be *symmetric*
129            uplos <- sapply(Tlst, slot, "uplo") ## either "U" or "L"
130            tLU <- table(uplos)# of length 1 or 2 ..
131            if(length(tLU) == 1) { ## all "U" or all "L"
132                useU <- uplos[1] == "U"
133            } else { ## length(tLU) == 2, counting "L" and "U"
134                useU <- diff(tLU) >= 0
135                if(useU && (hasL <- tLU[1] > 0))
136                    Tlst[hasL] <- lapply(Tlst[hasL], t)
137                else if(!useU && (hasU <- tLU[2] > 0))
138                    Tlst[hasU] <- lapply(Tlst[hasU], t)
139            }
140            if(use.n) { ## return nsparseMatrix :
141                r <- new("nsTMatrix")
142            } else {
143                r <- new("dsTMatrix")
144                r@x <- unlist(lapply(Tlst, slot, "x"))
145            }
146            r@uplo <- if(useU) "U" else "L"
147        }
148        else if(all(tri) && { ULs <- sapply(Tlst, slot, "uplo")##  "U" or "L"
149                              all(ULs[1L] == ULs[-1L]) } ## all upper or all lower
150           ){ ## *triangular* result
151    
152            if(use.n) { ## return nsparseMatrix :
153                r <- new("ntTMatrix")
154            } else {
155                r <- new("dtTMatrix")
156                r@x <- unlist(lapply(Tlst, slot, "x"))
157            }
158            r@uplo <- ULs[1L]
159        }
160        else {
161            if(any(sym))
162                Tlst[sym] <- lapply(Tlst[sym], as, "generalMatrix")
163            if(use.n) { ## return nsparseMatrix :
164                r <- new("ngTMatrix")
165            } else {
166                r <- new("dgTMatrix")
167                r@x <- unlist(lapply(Tlst, slot, "x"))
168            }
169        }
170        r@Dim <- c(i_off[nl+1], j_off[nl + 1])
171        r@i <- unlist(lapply(1:nl, function(k) Tlst[[k]]@i + i_off[k]))
172        r@j <- unlist(lapply(1:nl, function(k) Tlst[[k]]@j + j_off[k]))
173        r
174    }
175    
176    bdiag <- function(...) {
177        if((nA <- nargs()) == 0) return(new("dgCMatrix"))
178        if(nA == 1 && !is.list(...))
179            return(as(..., "CsparseMatrix"))
180        alis <- if(nA == 1 && is.list(..1)) ..1 else list(...)
181        if(length(alis) == 1)
182            return(as(alis[[1]], "CsparseMatrix"))
183    
184        ## else : two or more arguments
185        as(.bdiag(alis), "CsparseMatrix")
186  }  }
187    
188    
# Line 124  Line 233 
233  ## ddi*:  ## ddi*:
234  diag2tT <- function(from) .diag2tT(from, "U", "d")  diag2tT <- function(from) .diag2tT(from, "U", "d")
235  setAs("ddiMatrix", "triangularMatrix", diag2tT)  setAs("ddiMatrix", "triangularMatrix", diag2tT)
236  setAs("ddiMatrix", "sparseMatrix", diag2tT)  ##_no_longer_ setAs("ddiMatrix", "sparseMatrix", diag2tT)
237  ## needed too (otherwise <dense> -> Tsparse is taken):  ## needed too (otherwise <dense> -> Tsparse is taken):
238  setAs("ddiMatrix", "TsparseMatrix", diag2tT)  setAs("ddiMatrix", "TsparseMatrix", diag2tT)
239    setAs("ddiMatrix", "dsparseMatrix", diag2tT)
240  setAs("ddiMatrix", "CsparseMatrix",  setAs("ddiMatrix", "CsparseMatrix",
241        function(from) as(.diag2tT(from, "U", "d"), "CsparseMatrix"))        function(from) as(.diag2tT(from, "U", "d"), "CsparseMatrix"))
242  setAs("ddiMatrix", "symmetricMatrix",  setAs("ddiMatrix", "symmetricMatrix",
# Line 135  Line 245 
245  ## ldi*:  ## ldi*:
246  diag2tT <- function(from) .diag2tT(from, "U", "l")  diag2tT <- function(from) .diag2tT(from, "U", "l")
247  setAs("ldiMatrix", "triangularMatrix", diag2tT)  setAs("ldiMatrix", "triangularMatrix", diag2tT)
248  setAs("ldiMatrix", "sparseMatrix", diag2tT)  ##_no_longer_ setAs("ldiMatrix", "sparseMatrix", diag2tT)
249  ## needed too (otherwise <dense> -> Tsparse is taken):  ## needed too (otherwise <dense> -> Tsparse is taken):
250  setAs("ldiMatrix", "TsparseMatrix", diag2tT)  setAs("ldiMatrix", "TsparseMatrix", diag2tT)
251    setAs("ldiMatrix", "lsparseMatrix", diag2tT)
252  setAs("ldiMatrix", "CsparseMatrix",  setAs("ldiMatrix", "CsparseMatrix",
253        function(from) as(.diag2tT(from, "U", "l"), "CsparseMatrix"))        function(from) as(.diag2tT(from, "U", "l"), "CsparseMatrix"))
254  setAs("ldiMatrix", "symmetricMatrix",  setAs("ldiMatrix", "symmetricMatrix",
# Line 152  Line 263 
263                Dim = from@Dim, Dimnames = from@Dimnames)                Dim = from@Dim, Dimnames = from@Dimnames)
264        })        })
265    
266    setAs("diagonalMatrix", "nsparseMatrix", function(from) as(from, "nMatrix"))
267    
268  ## Cheap fast substitute for diag() which *does* preserve the mode of x :  ## Cheap fast substitute for diag() which *does* preserve the mode of x :
269  mkDiag <- function(x, n) {  mkDiag <- function(x, n) {
# Line 181  Line 293 
293  setAs("diagonalMatrix", "generalMatrix", # prefer sparse:  setAs("diagonalMatrix", "generalMatrix", # prefer sparse:
294        function(from) as(as(from, "CsparseMatrix"), "generalMatrix"))        function(from) as(as(from, "CsparseMatrix"), "generalMatrix"))
295    
296    setAs("diagonalMatrix", "denseMatrix",
297          function(from) as(as(from, "CsparseMatrix"), "denseMatrix"))
298    
299  .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
300    
301  .diag.2N <- function(m) {  .diag.2N <- function(m) {
# Line 201  Line 316 
316            d <- dim(from)            d <- dim(from)
317            if(d[1] != (n <- d[2])) stop("non-square matrix")            if(d[1] != (n <- d[2])) stop("non-square matrix")
318            if(any(from[row(from) != col(from)] != 0))            if(any(from[row(from) != col(from)] != 0))
319                stop("has non-zero off-diagonal entries")                stop("matrix with non-zero off-diagonals cannot be coerced to diagonalMatrix")
320            x <- diag(from)            x <- diag(from)
321            if(is.logical(x)) {            if(is.logical(x)) {
322                cl <- "ldiMatrix"                cl <- "ldiMatrix"
323                uni <- all(x)                uni <- allTrue(x) ## uni := {is it unit-diagonal ?}
324            } else {            } else {
325                cl <- "ddiMatrix"                cl <- "ddiMatrix"
326                uni <- all(x == 1)                uni <- allTrue(x == 1)
327                storage.mode(x) <- "double"                storage.mode(x) <- "double"
328            } ## TODO: complex            } ## TODO: complex
329            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 340 
340            x <- diag(from)            x <- diag(from)
341            if(is.logical(x)) {            if(is.logical(x)) {
342                cl <- "ldiMatrix"                cl <- "ldiMatrix"
343                uni <- all(x)                uni <- allTrue(x)
344            } else {            } else {
345                cl <- "ddiMatrix"                cl <- "ddiMatrix"
346                uni <- all(x == 1)                uni <- allTrue(x == 1)
347                storage.mode(x) <- "double"                storage.mode(x) <- "double"
348            }            } ## TODO: complex
349            new(cl, Dim = c(n,n), diag = if(uni) "U" else "N",            new(cl, Dim = c(n,n), diag = if(uni) "U" else "N",
350                x = if(uni) x[FALSE] else x)                x = if(uni) x[FALSE] else x)
351        })        })
# Line 240  Line 355 
355            function(x = 1, nrow, ncol) .diag.x(x))            function(x = 1, nrow, ncol) .diag.x(x))
356    
357  subDiag <- function(x, i, j, ..., drop) {  subDiag <- function(x, i, j, ..., drop) {
358      x <- as(x, "sparseMatrix")      x <- as(x, "CsparseMatrix") ## << was "TsparseMatrix" (Csparse is faster now)
359      x <- if(missing(i))      x <- if(missing(i))
360          x[, j, drop=drop]          x[, j, drop=drop]
361      else if(missing(j))      else if(missing(j))
362          x[i, , drop=drop]          if(nargs() == 4) x[i, , drop=drop] else x[i, drop=drop]
363      else      else
364          x[i,j, drop=drop]          x[i,j, drop=drop]
365      if(isS4(x) && isDiagonal(x)) as(x, "diagonalMatrix") else x      if(isS4(x) && isDiagonal(x)) as(x, "diagonalMatrix") else x
# Line 254  Line 369 
369                           j = "index", drop = "logical"), subDiag)                           j = "index", drop = "logical"), subDiag)
370  setMethod("[", signature(x = "diagonalMatrix", i = "index",  setMethod("[", signature(x = "diagonalMatrix", i = "index",
371                          j = "missing", drop = "logical"),                          j = "missing", drop = "logical"),
372            function(x, i, j, ..., drop) subDiag(x, i=i, drop=drop))            function(x, i, j, ..., drop) {
373                  na <- nargs()
374                  Matrix.msg("diag[i,m,l] : nargs()=", na, .M.level = 2)
375                  if(na == 4)
376                       subDiag(x, i=i, , drop=drop)
377                  else subDiag(x, i=i,   drop=drop)
378              })
379  setMethod("[", signature(x = "diagonalMatrix", i = "missing",  setMethod("[", signature(x = "diagonalMatrix", i = "missing",
380                           j = "index", drop = "logical"),                           j = "index", drop = "logical"),
381            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 385 
385  ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch  ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch
386  ## Only(?) current bug:  x[i] <- value  is wrong when  i is *vector*  ## Only(?) current bug:  x[i] <- value  is wrong when  i is *vector*
387  replDiag <- function(x, i, j, ..., value) {  replDiag <- function(x, i, j, ..., value) {
388      x <- as(x, "sparseMatrix")      x <- as(x, "TsparseMatrix")
389      if(missing(i))      if(missing(i))
390          x[, j] <- value          x[, j] <- value
391      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 430 
430                               x@x[ii] <- value                               x@x[ii] <- value
431                               x                               x
432                           } else { ## no longer diagonal, but remain sparse:                           } else { ## no longer diagonal, but remain sparse:
433                               x <- as(x, "sparseMatrix")                               x <- as(x, "TsparseMatrix")
434                               x[i] <- value                               x[i] <- value
435                               x                               x
436                           }                           }
# Line 325  Line 446 
446                                  j = "index", value = "replValue"),                                  j = "index", value = "replValue"),
447                   function(x,i,j, ..., value) replDiag(x, j=j, value=value))                   function(x,i,j, ..., value) replDiag(x, j=j, value=value))
448    
449    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",
450                                    value = "sparseMatrix"),
451                     function (x, i, j, ..., value)
452                     callGeneric(x=x, , j=j, value = as(value, "sparseVector")))
453    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "missing",
454                                    value = "sparseMatrix"),
455                     function (x, i, j, ..., value)
456                     callGeneric(x=x, i=i, , value = as(value, "sparseVector")))
457    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "index",
458                                    value = "sparseMatrix"),
459                     function (x, i, j, ..., value)
460                     callGeneric(x=x, i=i, j=j, value = as(value, "sparseVector")))
461    
462    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",
463                                    value = "sparseVector"),
464                     replDiag)
465    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "missing",
466                                    value = "sparseVector"),
467                     replDiag)
468    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "index",
469                                    value = "sparseVector"),
470                     replDiag)
471    
472    
473  setMethod("t", signature(x = "diagonalMatrix"),  setMethod("t", signature(x = "diagonalMatrix"),
474            function(x) { x@Dimnames <- x@Dimnames[2:1] ; x })            function(x) { x@Dimnames <- x@Dimnames[2:1] ; x })
# Line 371  Line 515 
515  ##       ---------------------  ##       ---------------------
516  ## Note that "ldi" logical are treated as numeric  ## Note that "ldi" logical are treated as numeric
517  diagdiagprod <- function(x, y) {  diagdiagprod <- function(x, y) {
518      n <- dimCheck(x,y)[1]      dimCheck(x,y)
519      if(x@diag != "U") {      if(x@diag != "U") {
520          if(y@diag != "U") {          if(y@diag != "U") {
521              nx <- x@x * y@x              nx <- x@x * y@x
# Line 403  Line 547 
547      dx <- dim(x)      dx <- dim(x)
548      dy <- dim(y)      dy <- dim(y)
549      if(dx[2] != dy[1]) stop("non-matching dimensions")      if(dx[2] != dy[1]) stop("non-matching dimensions")
     n <- dx[1]  
550      as(if(x@diag == "U") y else x@x * y, "Matrix")      as(if(x@diag == "U") y else x@x * y, "Matrix")
551  }  }
552  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),
# Line 472  Line 615 
615  ##        function(x, y = NULL) {  ##        function(x, y = NULL) {
616  ##           })  ##           })
617    
618  ## setMethod("tcrossprod", signature(x = "denseMatrix", y = "diagonalMatrix"),  Cspdiagprod <- function(x, y) {
619  ##        function(x, y = NULL) {      dx <- dim(x <- .Call(Csparse_diagU2N, x))
620  ##           })      dy <- dim(y)
621        if(dx[2] != dy[1]) stop("non-matching dimensions")
622        ind <- rep.int(seq_len(dx[2]), x@p[-1] - x@p[-dx[2]-1L])
623        if(y@diag == "N")
624            x@x <- x@x * y@x[ind]
625        if(is(x, "compMatrix") && length(xf <- x@factors)) {
626            ## instead of dropping all factors, be smart about some
627            ## TODO ......
628            x@factors <- list()
629        }
630        x
631    }
632    
633    diagCspprod <- function(x, y) {
634        dx <- dim(x)
635        dy <- dim(y <- .Call(Csparse_diagU2N, y))
636        if(dx[2] != dy[1]) stop("non-matching dimensions")
637        if(x@diag == "N")
638            y@x <- y@x * x@x[y@i + 1L]
639        if(is(y, "compMatrix") && length(yf <- y@factors)) {
640            ## instead of dropping all factors, be smart about some
641            ## TODO
642            keep <- character()
643            if(iLU <- names(yf) == "LU") {
644                ## TODO keep <- "LU"
645            }
646            y@factors <- yf[keep]
647        }
648        y
649    }
650    
651    setMethod("crossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
652              function(x, y = NULL) diagCspprod(x, y))
653    
654  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
655            function(x, y = NULL) crossprod(as(x, "sparseMatrix"), y))            function(x, y = NULL) diagCspprod(x, as(y, "CsparseMatrix")))
656    
657    ## Prefer calling diagCspprod to Cspdiagprod if going to transpose anyway
658    ##  x'y == (y'x)'
659    setMethod("crossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
660              function(x, y = NULL) t(diagCspprod(y, x)))
661    
662  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
663            function(x, y = NULL) crossprod(x, as(y, "sparseMatrix")))            function(x, y = NULL) t(diagCspprod(y, as(x, "Csparsematrix"))))
664    
665    setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
666              function(x, y = NULL) diagCspprod(x, t(y)))
667    
668  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
669            function(x, y = NULL) tcrossprod(as(x, "sparseMatrix"), y))            function(x, y = NULL) diagCspprod(x, t(as(y, "CsparseMatrix"))))
670    
671    setMethod("tcrossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
672              function(x, y = NULL) Cspdiagprod(x, y))
673    
674  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
675            function(x, y = NULL) tcrossprod(x, as(y, "sparseMatrix")))            function(x, y = NULL) Cspdiagprod(as(x, "CsparseMatrix"), y))
676    
677    setMethod("%*%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
678              function(x, y) diagCspprod(x, y))
679    
 ## FIXME?: In theory, this can be done *FASTER*, in some cases, via tapply1()  
680  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
681            function(x, y) as(x, "sparseMatrix") %*% y)            function(x, y) diagCspprod(as(x, "CsparseMatrix"), y))
682    
683  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
684            function(x, y) x %*% as(y, "sparseMatrix"))            function(x, y) Cspdiagprod(as(x, "CsparseMatrix"), y))
685  ## 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)  
686  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
687            function(x, y) x %*% as(y, "CsparseMatrix"))            function(x, y) Cspdiagprod(x, y))
688  ## NB: this is *not* needed for Tsparse & Rsparse  
689  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*
690  ##       do indeed work by going through sparse (and *not* ddense)!  ##       do indeed work by going through sparse (and *not* ddense)!
691    
   
   
692  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),
693            function(a, b, ...) {            function(a, b, ...) {
694                a@x <- 1/ a@x                a@x <- 1/ a@x
# Line 515  Line 697 
697            })            })
698    
699  solveDiag <- function(a, b, ...) {  solveDiag <- function(a, b, ...) {
700      if((n <- a@Dim[1]) != nrow(b))      if(a@Dim[1] != nrow(b))
701          stop("incompatible matrix dimensions")          stop("incompatible matrix dimensions")
702      ## trivially invert a 'in place' and multiply:      ## trivially invert a 'in place' and multiply:
703      a@x <- 1/ a@x      a@x <- 1/ a@x
# Line 556  Line 738 
738          .diag.2N(e1)          .diag.2N(e1)
739      }      }
740      else { ## result not diagonal, but at least symmetric:      else { ## result not diagonal, but at least symmetric:
741            ## e.g., m == m
742          isNum <- (is.numeric(r) || is.numeric(r00))          isNum <- (is.numeric(r) || is.numeric(r00))
743          isLog <- (is.logical(r) || is.logical(r00))          isLog <- (is.logical(r) || is.logical(r00))
744            Matrix.msg("exploding <diag> o <diag> into dense matrix", .M.level = 2)
         if(getOption("verbose"))  
             message("exploding  <diag>  o  <diag>  into dense matrix")  
745          d <- e1@Dim          d <- e1@Dim
746          n <- d[1]          n <- d[1]
747          stopifnot(length(r) == n)          stopifnot(length(r) == n)
# Line 588  Line 769 
769              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)
770  }  }
771    
772  ## FIXME:    diagonal  o  triangular  |-->  triangular  ## diagonal  o  triangular  |-->  triangular
773  ## -----     diagonal  o  symmetric   |-->  symmetric  ## diagonal  o  symmetric   |-->  symmetric
774  ##    {also when other is sparse: do these "here" --  ##    {also when other is sparse: do these "here" --
775  ##     before conversion to sparse, since that loses "diagonality"}  ##     before conversion to sparse, since that loses "diagonality"}
776    
# Line 619  Line 800 
800  ## ddi*: Arith: result numeric, potentially ddiMatrix  ## ddi*: Arith: result numeric, potentially ddiMatrix
801  setMethod("Arith", signature(e1 = "ddiMatrix", e2 = "numeric"),  setMethod("Arith", signature(e1 = "ddiMatrix", e2 = "numeric"),
802            function(e1,e2) {            function(e1,e2) {
803                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]; nsq <- n^2
804                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
805                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
806                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
# Line 637  Line 818 
818    
819  setMethod("Arith", signature(e1 = "numeric", e2 = "ddiMatrix"),  setMethod("Arith", signature(e1 = "numeric", e2 = "ddiMatrix"),
820            function(e1,e2) {            function(e1,e2) {
821                n <- e2@Dim[1]; nsq <- n*n                n <- e2@Dim[1]; nsq <- n^2
822                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
823                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
824                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
# Line 656  Line 837 
837  ## ldi* Arith --> result numeric, potentially ddiMatrix  ## ldi* Arith --> result numeric, potentially ddiMatrix
838  setMethod("Arith", signature(e1 = "ldiMatrix", e2 = "numeric"),  setMethod("Arith", signature(e1 = "ldiMatrix", e2 = "numeric"),
839            function(e1,e2) {            function(e1,e2) {
840                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]; nsq <- n^2
841                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
842                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
843                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
# Line 675  Line 856 
856    
857  setMethod("Arith", signature(e1 = "numeric", e2 = "ldiMatrix"),  setMethod("Arith", signature(e1 = "numeric", e2 = "ldiMatrix"),
858            function(e1,e2) {            function(e1,e2) {
859                n <- e2@Dim[1]; nsq <- n*n                n <- e2@Dim[1]; nsq <- n^2
860                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
861                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
862                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
# Line 695  Line 876 
876  ## ddi*: for "Ops" without Arith --> result logical, potentially ldi  ## ddi*: for "Ops" without Arith --> result logical, potentially ldi
877  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "numeric"),  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "numeric"),
878            function(e1,e2) {            function(e1,e2) {
879                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]; nsq <- n^2
880                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
881                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
882                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
# Line 714  Line 895 
895    
896  setMethod("Ops", signature(e1 = "numeric", e2 = "ddiMatrix"),  setMethod("Ops", signature(e1 = "numeric", e2 = "ddiMatrix"),
897            function(e1,e2) {            function(e1,e2) {
898                n <- e2@Dim[1]; nsq <- n*n                n <- e2@Dim[1]; nsq <- n^2
899                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
900                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
901                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
# Line 734  Line 915 
915  ## ldi*: for "Ops" without Arith --> result logical, potentially ldi  ## ldi*: for "Ops" without Arith --> result logical, potentially ldi
916  setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "numeric"),  setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "numeric"),
917            function(e1,e2) {            function(e1,e2) {
918                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]; nsq <- n^2
919                f0 <- callGeneric(FALSE, e2)                f0 <- callGeneric(FALSE, e2)
920                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
921                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
# Line 752  Line 933 
933    
934  setMethod("Ops", signature(e1 = "numeric", e2 = "ldiMatrix"),  setMethod("Ops", signature(e1 = "numeric", e2 = "ldiMatrix"),
935            function(e1,e2) {            function(e1,e2) {
936                n <- e2@Dim[1]; nsq <- n*n                n <- e2@Dim[1]; nsq <- n^2
937                f0 <- callGeneric(e1, FALSE)                f0 <- callGeneric(e1, FALSE)
938                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
939                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
# Line 774  Line 955 
955  for(other in c("ANY", "Matrix", "dMatrix")) {  for(other in c("ANY", "Matrix", "dMatrix")) {
956      ## ddi*:      ## ddi*:
957      setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other),      setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other),
958                function(e1,e2) callGeneric(diag2tT.u(e1,e2, "d"), e2))                function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2))
959      setMethod("Ops", signature(e1 = other, e2 = "ddiMatrix"),      setMethod("Ops", signature(e1 = other, e2 = "ddiMatrix"),
960                function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "d")))                function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "d")))
961      ## ldi*:      ## ldi*:
962      setMethod("Ops", signature(e1 = "ldiMatrix", e2 = other),      setMethod("Ops", signature(e1 = "ldiMatrix", e2 = other),
963                function(e1,e2) callGeneric(diag2tT.u(e1,e2, "l"), e2))                function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2))
964      setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),      setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),
965                function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "l")))                function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "l")))
966    }
967    
968    ## Direct subclasses of "denseMatrix": currently ddenseMatrix, ldense... :
969    if(FALSE)## too general, would contain  denseModelMatrix:
970    dense.subCl <- local({ dM.scl <- getClass("denseMatrix")@subclasses
971                           names(dM.scl)[sapply(dM.scl, slot, "distance") == 1] })
972    dense.subCl <- paste(c("d","l","n"), "denseMatrix", sep="")
973    for(DI in diCls) {
974        for(c2 in c(dense.subCl, "Matrix")) {
975            for(Fun in c("*", "^", "&")) {
976                setMethod(Fun, signature(e1 = DI, e2 = c2),
977                          function(e1,e2) callGeneric(e1, Diagonal(x = diag(e2))))
978                setMethod(Fun, signature(e1 = c2, e2 = DI),
979                          function(e1,e2) callGeneric(Diagonal(x = diag(e1)), e2))
980            }
981            ## NB: This arguably implicitly uses  0/0 :== 0  to keep diagonality
982            for(Fun in c("%%", "%/%", "/")) {
983                setMethod(Fun, signature(e1 = DI, e2 = c2),
984                          function(e1,e2) callGeneric(e1, Diagonal(x = diag(e2))))
985  }  }
 ## 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)  
986      }      }
987  }  }
988    
989    
990  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"
991  ### ----------  any, all: separately here  ### ----------   the last 4: separately here
992  for(cl in diCls) {  for(cl in diCls) {
993  setMethod("any", cl,  setMethod("any", cl,
994            function (x, ..., na.rm) {            function (x, ..., na.rm) {
995                if(any(x@Dim == 0)) FALSE                if(any(x@Dim == 0)) FALSE
996                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)
997            })            })
998  setMethod("all",  cl, function (x, ..., na.rm) any(x@Dim == 0))  setMethod("all",  cl, function (x, ..., na.rm) {
999  setMethod("prod", cl, function (x, ..., na.rm) as.numeric(any(x@Dim == 0)))      n <- x@Dim[1]
1000        if(n >= 2) FALSE
1001        else if(n == 0 || x@diag == "U") TRUE
1002        else all(x@x, ..., na.rm = na.rm)
1003    })
1004    setMethod("prod", cl, function (x, ..., na.rm) {
1005        n <- x@Dim[1]
1006        if(n >= 2) 0
1007        else if(n == 0 || x@diag == "U") 1
1008        else ## n == 1, diag = "N" :
1009            prod(x@x, ..., na.rm = na.rm)
1010    })
1011    
1012  setMethod("sum", cl,  setMethod("sum", cl,
1013            function(x, ..., na.rm) {            function(x, ..., na.rm) {
# Line 847  Line 1046 
1046      invisible(x)      invisible(x)
1047  }  }
1048    
1049    ## somewhat consistent with "print" for sparseMatrix :
1050    setMethod("print", signature(x = "diagonalMatrix"), prDiag)
1051    
1052  setMethod("show", signature(object = "diagonalMatrix"),  setMethod("show", signature(object = "diagonalMatrix"),
1053            function(object) {            function(object) {
1054                d <- dim(object)                d <- dim(object)
1055                cl <- class(object)                cl <- class(object)
1056                cat(sprintf('%d x %d diagonal matrix of class "%s"\n',                cat(sprintf('%d x %d diagonal matrix of class "%s"',
1057                            d[1], d[2], cl))                            d[1], d[2], cl))
1058                  if(d[1] < 50) {
1059                      cat("\n")
1060                prDiag(object)                prDiag(object)
1061                  } else {
1062                      cat(", with diagonal entries\n")
1063                      show(diag(object))
1064                      invisible(object)
1065                  }
1066            })            })

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

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