SCM

SCM Repository

[matrix] Diff of /pkg/Matrix/R/diagMatrix.R
ViewVC logotype

Diff of /pkg/Matrix/R/diagMatrix.R

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 2197, Mon Jun 2 14:34:42 2008 UTC revision 2544, Wed Jun 2 09:29:01 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        x
627    }
628    
629    diagCspprod <- function(x, y) {
630        dx <- dim(x)
631        dy <- dim(y <- .Call(Csparse_diagU2N, y))
632        if(dx[2] != dy[1]) stop("non-matching dimensions")
633        if(x@diag == "N")
634            y@x <- y@x * x@x[y@i + 1L]
635        y
636    }
637    
638    setMethod("crossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
639              function(x, y = NULL) diagCspprod(x, y))
640    
641  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
642            function(x, y = NULL) crossprod(as(x, "sparseMatrix"), y))            function(x, y = NULL) diagCspprod(x, as(y, "CsparseMatrix")))
643    
644    ## Prefer calling diagCspprod to Cspdiagprod if going to transpose anyway
645    ##  x'y == (y'x)'
646    setMethod("crossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
647              function(x, y = NULL) t(diagCspprod(y, x)))
648    
649  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
650            function(x, y = NULL) crossprod(x, as(y, "sparseMatrix")))            function(x, y = NULL) t(diagCspprod(y, as(x, "Csparsematrix"))))
651    
652    setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
653              function(x, y = NULL) diagCspprod(x, t(y)))
654    
655  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
656            function(x, y = NULL) tcrossprod(as(x, "sparseMatrix"), y))            function(x, y = NULL) diagCspprod(x, t(as(y, "CsparseMatrix"))))
657    
658    setMethod("tcrossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
659              function(x, y = NULL) Cspdiagprod(x, y))
660    
661  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
662            function(x, y = NULL) tcrossprod(x, as(y, "sparseMatrix")))            function(x, y = NULL) Cspdiagprod(as(x, "CsparseMatrix"), y))
663    
664    setMethod("%*%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
665              function(x, y) diagCspprod(x, y))
666    
 ## FIXME?: In theory, this can be done *FASTER*, in some cases, via tapply1()  
667  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
668            function(x, y) as(x, "sparseMatrix") %*% y)            function(x, y) diagCspprod(as(x, "CsparseMatrix"), y))
669    
670  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
671            function(x, y) x %*% as(y, "sparseMatrix"))            function(x, y) Cspdiagprod(as(x, "CsparseMatrix"), y))
672  ## 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)  
673  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
674            function(x, y) x %*% as(y, "CsparseMatrix"))            function(x, y) Cspdiagprod(x, y))
675  ## NB: this is *not* needed for Tsparse & Rsparse  
676  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*
677  ##       do indeed work by going through sparse (and *not* ddense)!  ##       do indeed work by going through sparse (and *not* ddense)!
678    
   
   
679  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),
680            function(a, b, ...) {            function(a, b, ...) {
681                a@x <- 1/ a@x                a@x <- 1/ a@x
# Line 556  Line 725 
725          .diag.2N(e1)          .diag.2N(e1)
726      }      }
727      else { ## result not diagonal, but at least symmetric:      else { ## result not diagonal, but at least symmetric:
728            ## e.g., m == m
729          isNum <- (is.numeric(r) || is.numeric(r00))          isNum <- (is.numeric(r) || is.numeric(r00))
730          isLog <- (is.logical(r) || is.logical(r00))          isLog <- (is.logical(r) || is.logical(r00))
731            Matrix.msg("exploding <diag> o <diag> into dense matrix", .M.level = 2)
         if(getOption("verbose"))  
             message("exploding  <diag>  o  <diag>  into dense matrix")  
732          d <- e1@Dim          d <- e1@Dim
733          n <- d[1]          n <- d[1]
734          stopifnot(length(r) == n)          stopifnot(length(r) == n)
# Line 588  Line 756 
756              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)
757  }  }
758    
759  ## FIXME:    diagonal  o  triangular  |-->  triangular  ## diagonal  o  triangular  |-->  triangular
760  ## -----     diagonal  o  symmetric   |-->  symmetric  ## diagonal  o  symmetric   |-->  symmetric
761  ##    {also when other is sparse: do these "here" --  ##    {also when other is sparse: do these "here" --
762  ##     before conversion to sparse, since that loses "diagonality"}  ##     before conversion to sparse, since that loses "diagonality"}
763    
# Line 619  Line 787 
787  ## ddi*: Arith: result numeric, potentially ddiMatrix  ## ddi*: Arith: result numeric, potentially ddiMatrix
788  setMethod("Arith", signature(e1 = "ddiMatrix", e2 = "numeric"),  setMethod("Arith", signature(e1 = "ddiMatrix", e2 = "numeric"),
789            function(e1,e2) {            function(e1,e2) {
790                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]; nsq <- n^2
791                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
792                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
793                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
# Line 637  Line 805 
805    
806  setMethod("Arith", signature(e1 = "numeric", e2 = "ddiMatrix"),  setMethod("Arith", signature(e1 = "numeric", e2 = "ddiMatrix"),
807            function(e1,e2) {            function(e1,e2) {
808                n <- e2@Dim[1]; nsq <- n*n                n <- e2@Dim[1]; nsq <- n^2
809                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
810                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
811                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
# Line 656  Line 824 
824  ## ldi* Arith --> result numeric, potentially ddiMatrix  ## ldi* Arith --> result numeric, potentially ddiMatrix
825  setMethod("Arith", signature(e1 = "ldiMatrix", e2 = "numeric"),  setMethod("Arith", signature(e1 = "ldiMatrix", e2 = "numeric"),
826            function(e1,e2) {            function(e1,e2) {
827                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]; nsq <- n^2
828                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
829                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
830                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
# Line 675  Line 843 
843    
844  setMethod("Arith", signature(e1 = "numeric", e2 = "ldiMatrix"),  setMethod("Arith", signature(e1 = "numeric", e2 = "ldiMatrix"),
845            function(e1,e2) {            function(e1,e2) {
846                n <- e2@Dim[1]; nsq <- n*n                n <- e2@Dim[1]; nsq <- n^2
847                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
848                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
849                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
# Line 695  Line 863 
863  ## ddi*: for "Ops" without Arith --> result logical, potentially ldi  ## ddi*: for "Ops" without Arith --> result logical, potentially ldi
864  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "numeric"),  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "numeric"),
865            function(e1,e2) {            function(e1,e2) {
866                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]; nsq <- n^2
867                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
868                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
869                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
# Line 714  Line 882 
882    
883  setMethod("Ops", signature(e1 = "numeric", e2 = "ddiMatrix"),  setMethod("Ops", signature(e1 = "numeric", e2 = "ddiMatrix"),
884            function(e1,e2) {            function(e1,e2) {
885                n <- e2@Dim[1]; nsq <- n*n                n <- e2@Dim[1]; nsq <- n^2
886                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
887                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
888                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
# Line 734  Line 902 
902  ## ldi*: for "Ops" without Arith --> result logical, potentially ldi  ## ldi*: for "Ops" without Arith --> result logical, potentially ldi
903  setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "numeric"),  setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "numeric"),
904            function(e1,e2) {            function(e1,e2) {
905                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]; nsq <- n^2
906                f0 <- callGeneric(FALSE, e2)                f0 <- callGeneric(FALSE, e2)
907                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
908                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
# Line 752  Line 920 
920    
921  setMethod("Ops", signature(e1 = "numeric", e2 = "ldiMatrix"),  setMethod("Ops", signature(e1 = "numeric", e2 = "ldiMatrix"),
922            function(e1,e2) {            function(e1,e2) {
923                n <- e2@Dim[1]; nsq <- n*n                n <- e2@Dim[1]; nsq <- n^2
924                f0 <- callGeneric(e1, FALSE)                f0 <- callGeneric(e1, FALSE)
925                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
926                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
# Line 774  Line 942 
942  for(other in c("ANY", "Matrix", "dMatrix")) {  for(other in c("ANY", "Matrix", "dMatrix")) {
943      ## ddi*:      ## ddi*:
944      setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other),      setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other),
945                function(e1,e2) callGeneric(diag2tT.u(e1,e2, "d"), e2))                function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2))
946      setMethod("Ops", signature(e1 = other, e2 = "ddiMatrix"),      setMethod("Ops", signature(e1 = other, e2 = "ddiMatrix"),
947                function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "d")))                function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "d")))
948      ## ldi*:      ## ldi*:
949      setMethod("Ops", signature(e1 = "ldiMatrix", e2 = other),      setMethod("Ops", signature(e1 = "ldiMatrix", e2 = other),
950                function(e1,e2) callGeneric(diag2tT.u(e1,e2, "l"), e2))                function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2))
951      setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),      setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),
952                function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "l")))                function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "l")))
953    }
954    
955    ## Direct subclasses of "denseMatrix": currently ddenseMatrix, ldense... :
956    if(FALSE)## too general, would contain  denseModelMatrix:
957    dense.subCl <- local({ dM.scl <- getClass("denseMatrix")@subclasses
958                           names(dM.scl)[sapply(dM.scl, slot, "distance") == 1] })
959    dense.subCl <- paste(c("d","l","n"), "denseMatrix", sep="")
960    for(DI in diCls) {
961        for(c2 in c(dense.subCl, "Matrix")) {
962            for(Fun in c("*", "^", "&")) {
963                setMethod(Fun, signature(e1 = DI, e2 = c2),
964                          function(e1,e2) callGeneric(e1, Diagonal(x = diag(e2))))
965                setMethod(Fun, signature(e1 = c2, e2 = DI),
966                          function(e1,e2) callGeneric(Diagonal(x = diag(e1)), e2))
967            }
968            ## NB: This arguably implicitly uses  0/0 :== 0  to keep diagonality
969            for(Fun in c("%%", "%/%", "/")) {
970                setMethod(Fun, signature(e1 = DI, e2 = c2),
971                          function(e1,e2) callGeneric(e1, Diagonal(x = diag(e2))))
972  }  }
 ## 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)  
973      }      }
974  }  }
975    
976    
977  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"
978  ### ----------  any, all: separately here  ### ----------   the last 4: separately here
979  for(cl in diCls) {  for(cl in diCls) {
980  setMethod("any", cl,  setMethod("any", cl,
981            function (x, ..., na.rm) {            function (x, ..., na.rm) {
982                if(any(x@Dim == 0)) FALSE                if(any(x@Dim == 0)) FALSE
983                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)
984            })            })
985  setMethod("all",  cl, function (x, ..., na.rm) any(x@Dim == 0))  setMethod("all",  cl, function (x, ..., na.rm) {
986  setMethod("prod", cl, function (x, ..., na.rm) as.numeric(any(x@Dim == 0)))      n <- x@Dim[1]
987        if(n >= 2) FALSE
988        else if(n == 0 || x@diag == "U") TRUE
989        else all(x@x, ..., na.rm = na.rm)
990    })
991    setMethod("prod", cl, function (x, ..., na.rm) {
992        n <- x@Dim[1]
993        if(n >= 2) 0
994        else if(n == 0 || x@diag == "U") 1
995        else ## n == 1, diag = "N" :
996            prod(x@x, ..., na.rm = na.rm)
997    })
998    
999  setMethod("sum", cl,  setMethod("sum", cl,
1000            function(x, ..., na.rm) {            function(x, ..., na.rm) {
# Line 847  Line 1033 
1033      invisible(x)      invisible(x)
1034  }  }
1035    
1036    ## somewhat consistent with "print" for sparseMatrix :
1037    setMethod("print", signature(x = "diagonalMatrix"), prDiag)
1038    
1039  setMethod("show", signature(object = "diagonalMatrix"),  setMethod("show", signature(object = "diagonalMatrix"),
1040            function(object) {            function(object) {
1041                d <- dim(object)                d <- dim(object)
1042                cl <- class(object)                cl <- class(object)
1043                cat(sprintf('%d x %d diagonal matrix of class "%s"\n',                cat(sprintf('%d x %d diagonal matrix of class "%s"',
1044                            d[1], d[2], cl))                            d[1], d[2], cl))
1045                  if(d[1] < 50) {
1046                      cat("\n")
1047                prDiag(object)                prDiag(object)
1048                  } else {
1049                      cat(", with diagonal entries\n")
1050                      show(diag(object))
1051                      invisible(object)
1052                  }
1053            })            })

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

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