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 2608, Tue Aug 10 08:23:48 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 472  Line 616 
616  ##        function(x, y = NULL) {  ##        function(x, y = NULL) {
617  ##           })  ##           })
618    
619  ## setMethod("tcrossprod", signature(x = "denseMatrix", y = "diagonalMatrix"),  Cspdiagprod <- function(x, y) {
620  ##        function(x, y = NULL) {      dx <- dim(x <- .Call(Csparse_diagU2N, x))
621  ##           })      dy <- dim(y)
622        if(dx[2] != dy[1]) stop("non-matching dimensions")
623        ind <- rep.int(seq_len(dx[2]), x@p[-1] - x@p[-dx[2]-1L])
624        if(y@diag == "N")
625            x@x <- x@x * y@x[ind]
626        if(is(x, "compMatrix") && length(xf <- x@factors)) {
627            ## instead of dropping all factors, be smart about some
628            ## TODO ......
629            x@factors <- list()
630        }
631        x
632    }
633    
634    diagCspprod <- function(x, y) {
635        dx <- dim(x)
636        dy <- dim(y <- .Call(Csparse_diagU2N, y))
637        if(dx[2] != dy[1]) stop("non-matching dimensions")
638        if(x@diag == "N")
639            y@x <- y@x * x@x[y@i + 1L]
640        if(is(y, "compMatrix") && length(yf <- y@factors)) {
641            ## instead of dropping all factors, be smart about some
642            ## TODO
643            keep <- character()
644            if(iLU <- names(yf) == "LU") {
645                ## TODO keep <- "LU"
646            }
647            y@factors <- yf[keep]
648        }
649        y
650    }
651    
652    setMethod("crossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
653              function(x, y = NULL) diagCspprod(x, y))
654    
655  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
656            function(x, y = NULL) crossprod(as(x, "sparseMatrix"), y))            function(x, y = NULL) diagCspprod(x, as(y, "CsparseMatrix")))
657    
658    ## Prefer calling diagCspprod to Cspdiagprod if going to transpose anyway
659    ##  x'y == (y'x)'
660    setMethod("crossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
661              function(x, y = NULL) t(diagCspprod(y, x)))
662    
663  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
664            function(x, y = NULL) crossprod(x, as(y, "sparseMatrix")))            function(x, y = NULL) t(diagCspprod(y, as(x, "Csparsematrix"))))
665    
666    setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
667              function(x, y = NULL) diagCspprod(x, t(y)))
668    
669  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
670            function(x, y = NULL) tcrossprod(as(x, "sparseMatrix"), y))            function(x, y = NULL) diagCspprod(x, t(as(y, "CsparseMatrix"))))
671    
672    setMethod("tcrossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
673              function(x, y = NULL) Cspdiagprod(x, y))
674    
675  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
676            function(x, y = NULL) tcrossprod(x, as(y, "sparseMatrix")))            function(x, y = NULL) Cspdiagprod(as(x, "CsparseMatrix"), y))
677    
678    setMethod("%*%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
679              function(x, y) diagCspprod(x, y))
680    
 ## FIXME?: In theory, this can be done *FASTER*, in some cases, via tapply1()  
681  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
682            function(x, y) as(x, "sparseMatrix") %*% y)            function(x, y) diagCspprod(as(x, "CsparseMatrix"), y))
683    
684  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
685            function(x, y) x %*% as(y, "sparseMatrix"))            function(x, y) Cspdiagprod(as(x, "CsparseMatrix"), y))
686  ## 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)  
687  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
688            function(x, y) x %*% as(y, "CsparseMatrix"))            function(x, y) Cspdiagprod(x, y))
689  ## NB: this is *not* needed for Tsparse & Rsparse  
690  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*
691  ##       do indeed work by going through sparse (and *not* ddense)!  ##       do indeed work by going through sparse (and *not* ddense)!
692    
   
   
693  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),
694            function(a, b, ...) {            function(a, b, ...) {
695                a@x <- 1/ a@x                a@x <- 1/ a@x
# Line 556  Line 739 
739          .diag.2N(e1)          .diag.2N(e1)
740      }      }
741      else { ## result not diagonal, but at least symmetric:      else { ## result not diagonal, but at least symmetric:
742            ## e.g., m == m
743          isNum <- (is.numeric(r) || is.numeric(r00))          isNum <- (is.numeric(r) || is.numeric(r00))
744          isLog <- (is.logical(r) || is.logical(r00))          isLog <- (is.logical(r) || is.logical(r00))
745            Matrix.msg("exploding <diag> o <diag> into dense matrix", .M.level = 2)
         if(getOption("verbose"))  
             message("exploding  <diag>  o  <diag>  into dense matrix")  
746          d <- e1@Dim          d <- e1@Dim
747          n <- d[1]          n <- d[1]
748          stopifnot(length(r) == n)          stopifnot(length(r) == n)
# Line 588  Line 770 
770              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)
771  }  }
772    
773  ## FIXME:    diagonal  o  triangular  |-->  triangular  ## diagonal  o  triangular  |-->  triangular
774  ## -----     diagonal  o  symmetric   |-->  symmetric  ## diagonal  o  symmetric   |-->  symmetric
775  ##    {also when other is sparse: do these "here" --  ##    {also when other is sparse: do these "here" --
776  ##     before conversion to sparse, since that loses "diagonality"}  ##     before conversion to sparse, since that loses "diagonality"}
777    
# Line 619  Line 801 
801  ## ddi*: Arith: result numeric, potentially ddiMatrix  ## ddi*: Arith: result numeric, potentially ddiMatrix
802  setMethod("Arith", signature(e1 = "ddiMatrix", e2 = "numeric"),  setMethod("Arith", signature(e1 = "ddiMatrix", e2 = "numeric"),
803            function(e1,e2) {            function(e1,e2) {
804                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]; nsq <- n^2
805                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
806                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
807                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
# Line 637  Line 819 
819    
820  setMethod("Arith", signature(e1 = "numeric", e2 = "ddiMatrix"),  setMethod("Arith", signature(e1 = "numeric", e2 = "ddiMatrix"),
821            function(e1,e2) {            function(e1,e2) {
822                n <- e2@Dim[1]; nsq <- n*n                n <- e2@Dim[1]; nsq <- n^2
823                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
824                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
825                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
# Line 656  Line 838 
838  ## ldi* Arith --> result numeric, potentially ddiMatrix  ## ldi* Arith --> result numeric, potentially ddiMatrix
839  setMethod("Arith", signature(e1 = "ldiMatrix", e2 = "numeric"),  setMethod("Arith", signature(e1 = "ldiMatrix", e2 = "numeric"),
840            function(e1,e2) {            function(e1,e2) {
841                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]; nsq <- n^2
842                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
843                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
844                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
# Line 675  Line 857 
857    
858  setMethod("Arith", signature(e1 = "numeric", e2 = "ldiMatrix"),  setMethod("Arith", signature(e1 = "numeric", e2 = "ldiMatrix"),
859            function(e1,e2) {            function(e1,e2) {
860                n <- e2@Dim[1]; nsq <- n*n                n <- e2@Dim[1]; nsq <- n^2
861                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
862                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
863                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
# Line 695  Line 877 
877  ## ddi*: for "Ops" without Arith --> result logical, potentially ldi  ## ddi*: for "Ops" without Arith --> result logical, potentially ldi
878  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "numeric"),  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "numeric"),
879            function(e1,e2) {            function(e1,e2) {
880                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]; nsq <- n^2
881                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
882                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
883                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
# Line 714  Line 896 
896    
897  setMethod("Ops", signature(e1 = "numeric", e2 = "ddiMatrix"),  setMethod("Ops", signature(e1 = "numeric", e2 = "ddiMatrix"),
898            function(e1,e2) {            function(e1,e2) {
899                n <- e2@Dim[1]; nsq <- n*n                n <- e2@Dim[1]; nsq <- n^2
900                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
901                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
902                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
# Line 734  Line 916 
916  ## ldi*: for "Ops" without Arith --> result logical, potentially ldi  ## ldi*: for "Ops" without Arith --> result logical, potentially ldi
917  setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "numeric"),  setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "numeric"),
918            function(e1,e2) {            function(e1,e2) {
919                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]; nsq <- n^2
920                f0 <- callGeneric(FALSE, e2)                f0 <- callGeneric(FALSE, e2)
921                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
922                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
# Line 752  Line 934 
934    
935  setMethod("Ops", signature(e1 = "numeric", e2 = "ldiMatrix"),  setMethod("Ops", signature(e1 = "numeric", e2 = "ldiMatrix"),
936            function(e1,e2) {            function(e1,e2) {
937                n <- e2@Dim[1]; nsq <- n*n                n <- e2@Dim[1]; nsq <- n^2
938                f0 <- callGeneric(e1, FALSE)                f0 <- callGeneric(e1, FALSE)
939                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
940                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
# Line 774  Line 956 
956  for(other in c("ANY", "Matrix", "dMatrix")) {  for(other in c("ANY", "Matrix", "dMatrix")) {
957      ## ddi*:      ## ddi*:
958      setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other),      setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other),
959                function(e1,e2) callGeneric(diag2tT.u(e1,e2, "d"), e2))                function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2))
960      setMethod("Ops", signature(e1 = other, e2 = "ddiMatrix"),      setMethod("Ops", signature(e1 = other, e2 = "ddiMatrix"),
961                function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "d")))                function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "d")))
962      ## ldi*:      ## ldi*:
963      setMethod("Ops", signature(e1 = "ldiMatrix", e2 = other),      setMethod("Ops", signature(e1 = "ldiMatrix", e2 = other),
964                function(e1,e2) callGeneric(diag2tT.u(e1,e2, "l"), e2))                function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2))
965      setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),      setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),
966                function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "l")))                function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "l")))
967    }
968    
969    ## Direct subclasses of "denseMatrix": currently ddenseMatrix, ldense... :
970    if(FALSE)## too general, would contain  denseModelMatrix:
971    dense.subCl <- local({ dM.scl <- getClass("denseMatrix")@subclasses
972                           names(dM.scl)[sapply(dM.scl, slot, "distance") == 1] })
973    dense.subCl <- paste(c("d","l","n"), "denseMatrix", sep="")
974    for(DI in diCls) {
975        for(c2 in c(dense.subCl, "Matrix")) {
976            for(Fun in c("*", "^", "&")) {
977                setMethod(Fun, signature(e1 = DI, e2 = c2),
978                          function(e1,e2) callGeneric(e1, Diagonal(x = diag(e2))))
979                setMethod(Fun, signature(e1 = c2, e2 = DI),
980                          function(e1,e2) callGeneric(Diagonal(x = diag(e1)), e2))
981            }
982            ## NB: This arguably implicitly uses  0/0 :== 0  to keep diagonality
983            for(Fun in c("%%", "%/%", "/")) {
984                setMethod(Fun, signature(e1 = DI, e2 = c2),
985                          function(e1,e2) callGeneric(e1, Diagonal(x = diag(e2))))
986  }  }
 ## 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)  
987      }      }
988  }  }
989    
990    
991  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"
992  ### ----------  any, all: separately here  ### ----------   the last 4: separately here
993  for(cl in diCls) {  for(cl in diCls) {
994  setMethod("any", cl,  setMethod("any", cl,
995            function (x, ..., na.rm) {            function (x, ..., na.rm) {
996                if(any(x@Dim == 0)) FALSE                if(any(x@Dim == 0)) FALSE
997                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)
998            })            })
999  setMethod("all",  cl, function (x, ..., na.rm) any(x@Dim == 0))  setMethod("all",  cl, function (x, ..., na.rm) {
1000  setMethod("prod", cl, function (x, ..., na.rm) as.numeric(any(x@Dim == 0)))      n <- x@Dim[1]
1001        if(n >= 2) FALSE
1002        else if(n == 0 || x@diag == "U") TRUE
1003        else all(x@x, ..., na.rm = na.rm)
1004    })
1005    setMethod("prod", cl, function (x, ..., na.rm) {
1006        n <- x@Dim[1]
1007        if(n >= 2) 0
1008        else if(n == 0 || x@diag == "U") 1
1009        else ## n == 1, diag = "N" :
1010            prod(x@x, ..., na.rm = na.rm)
1011    })
1012    
1013  setMethod("sum", cl,  setMethod("sum", cl,
1014            function(x, ..., na.rm) {            function(x, ..., na.rm) {
# Line 847  Line 1047 
1047      invisible(x)      invisible(x)
1048  }  }
1049    
1050    ## somewhat consistent with "print" for sparseMatrix :
1051    setMethod("print", signature(x = "diagonalMatrix"), prDiag)
1052    
1053  setMethod("show", signature(object = "diagonalMatrix"),  setMethod("show", signature(object = "diagonalMatrix"),
1054            function(object) {            function(object) {
1055                d <- dim(object)                d <- dim(object)
1056                cl <- class(object)                cl <- class(object)
1057                cat(sprintf('%d x %d diagonal matrix of class "%s"\n',                cat(sprintf('%d x %d diagonal matrix of class "%s"',
1058                            d[1], d[2], cl))                            d[1], d[2], cl))
1059                  if(d[1] < 50) {
1060                      cat("\n")
1061                prDiag(object)                prDiag(object)
1062                  } else {
1063                      cat(", with diagonal entries\n")
1064                      show(diag(object))
1065                      invisible(object)
1066                  }
1067            })            })

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

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