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 2820, Mon Aug 20 14:06:23 2012 UTC
# Line 5  Line 5 
5  ##          but *not* diag() extractor!  ##          but *not* diag() extractor!
6  Diagonal <- function(n, x = NULL)  Diagonal <- function(n, x = NULL)
7  {  {
8      ## Allow  Diagonal(4)  and  Diagonal(x=1:5)      ## Allow  Diagonal(4), Diagonal(x=1:5), and  Diagonal(4, TRUE)
9      if(missing(n))      if(missing(n))
10          n <- length(x)          n <- length(x)
11      else {      else {
# Line 17  Line 17 
17          new("ddiMatrix", Dim = c(n,n), diag = "U")          new("ddiMatrix", Dim = c(n,n), diag = "U")
18      else {      else {
19          lx <- length(x)          lx <- length(x)
20          stopifnot(lx == 1 || lx == n) # but keep 'x' short for now          lx.1 <- lx == 1L
21            stopifnot(lx.1 || lx == n) # but keep 'x' short for now
22          if(is.logical(x))          if(is.logical(x))
23              cl <- "ldiMatrix"              cl <- "ldiMatrix"
24          else if(is.numeric(x)) {          else if(is.numeric(x)) {
# Line 27  Line 28 
28          else if(is.complex(x)) {          else if(is.complex(x)) {
29              cl <- "zdiMatrix"  # will not yet work              cl <- "zdiMatrix"  # will not yet work
30          } else stop("'x' has invalid data type")          } else stop("'x' has invalid data type")
31            if(lx.1 && !is.na(x) && x == 1) # cheap check for uni-diagonal..
32                new(cl, Dim = c(n,n), diag = "U")
33            else
34          new(cl, Dim = c(n,n), diag = "N",          new(cl, Dim = c(n,n), diag = "N",
35              x = if(lx == 1) rep.int(x,n) else x)                  x = if(lx.1) rep.int(x,n) else x)
36      }      }
37  }  }
38    
39  ## Pkg 'spdep' had (relatively slow) versions of this as_dsCMatrix_I()  .sparseDiagonal <- function(n, x = rep.int(1,m), uplo = "U",
40  .symDiagonal <- function(n, x = rep.int(1,n), uplo = "U") {                              shape = if(missing(cols)) "t" else "g",
41                                kind, cols = if(n) 0:(n - 1L) else integer(0))
42    {
43      stopifnot(n == (n. <- as.integer(n)), (n <- n.) >= 0)      stopifnot(n == (n. <- as.integer(n)), (n <- n.) >= 0)
44      if((lx <- length(x)) == 1) x <- rep.int(x, n)      if(!missing(cols))
45      else if(lx != n) stop("length(x) must be 1 or n")          stopifnot(0 <= (cols <- as.integer(cols)), cols < n)
46      cls <-      m <- length(cols)
47          if(is.double(x)) "dsCMatrix"      if(missing(kind))
48          else if(is.logical(x)) "lsCMatrix"          kind <-
49                if(is.double(x)) "d"
50                else if(is.logical(x)) "l"
51          else { ## for now          else { ## for now
52              storage.mode(x) <- "double"              storage.mode(x) <- "double"
53              "dsCMatrix"                  "d"
54          }          }
55      new(cls, Dim = c(n,n), x = x, uplo = uplo,      else stopifnot(any(kind == c("d","l","n")))
56          i = if(n) 0:(n - 1L) else integer(0), p = 0:n)      if(kind != "n") {
57            if((lx <- length(x)) == 1) x <- rep.int(x, m)
58            else if(lx != m) stop("length(x) must be either 1 or #{cols}")
59        }
60        stopifnot(is.character(shape), nchar(shape) == 1,
61                  any(shape == c("t","s","g"))) # triangular / symmetric / general
62        if(kind == "n") {
63            if(shape == "g")
64                new("ngCMatrix", Dim = c(n,m), i = cols, p = 0:m)
65            else new(paste0("n", shape, "CMatrix"), Dim = c(n,m), uplo = uplo,
66                     i = cols, p = 0:m)
67        }
68        ## kind != "n" -- have x slot :
69        else if(shape == "g")
70            new(paste0(kind, "gCMatrix"), Dim = c(n,m),
71                x = x, i = cols, p = 0:m)
72        else new(paste0(kind, shape, "CMatrix"), Dim = c(n,m), uplo = uplo,
73                 x = x, i = cols, p = 0:m)
74  }  }
75    
76  ### This is modified from a post of Bert Gunter to R-help on  1 Sep 2005.  ## Pkg 'spdep' had (relatively slow) versions of this as_dsCMatrix_I()
77  ### Bert's code built on a post by Andy Liaw who most probably was influenced  .symDiagonal <- function(n, x = rep.int(1,n), uplo = "U")
78  ### by earlier posts, notably one by Scott Chasalow on S-news, 16 Jan 2002      .sparseDiagonal(n, x, uplo, shape = "s")
 ### who posted his bdiag() function written in December 1995.  
79    
80  bdiag <- function(...) {  # instead of   diagU2N(as(Diagonal(n), "CsparseMatrix")), diag = "N" in any case:
81      if(nargs() == 0) return(new("dgCMatrix"))  .trDiagonal <- function(n, x = rep.int(1,n), uplo = "U")
82      ## else :      .sparseDiagonal(n, x, uplo, shape = "t")
83      mlist <- if (nargs() == 1) as.list(...) else list(...)  
84      dims <- sapply(mlist, dim)  
85    ## This is modified from a post of Bert Gunter to R-help on  1 Sep 2005.
86    ## Bert's code built on a post by Andy Liaw who most probably was influenced
87    ## by earlier posts, notably one by Scott Chasalow on S-news, 16 Jan 2002
88    ## who posted his bdiag() function written in December 1995.
89    if(FALSE)##--- no longer used:
90    .bdiag <- function(lst) {
91        ## block-diagonal matrix [a dgTMatrix] from list of matrices
92        stopifnot(is.list(lst), length(lst) >= 1)
93        dims <- vapply(lst, dim, 1L, USE.NAMES=FALSE)
94      ## make sure we had all matrices:      ## make sure we had all matrices:
95      if(!(is.matrix(dims) && nrow(dims) == 2))      if(!(is.matrix(dims) && nrow(dims) == 2))
96          stop("some arguments are not matrices")          stop("some arguments are not matrices")
97      csdim <- rbind(rep.int(0L, 2),      csdim <- rbind(rep.int(0L, 2),
98                     apply(sapply(mlist, dim), 1, cumsum))                     apply(dims, 1, cumsum))
99      ret <- new("dgTMatrix", Dim = as.integer(csdim[nrow(csdim),]))      r <- new("dgTMatrix")
100        r@Dim <- as.integer(csdim[nrow(csdim),])
101      add1 <- matrix(1:0, 2,2)      add1 <- matrix(1:0, 2,2)
102      for(i in seq_along(mlist)) {      for(i in seq_along(lst)) {
103          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])
104          if(is.null(dim(indx))) ## non-square matrix          if(is.null(dim(indx))) ## non-square matrix
105              ret[indx[[1]],indx[[2]]] <- mlist[[i]]              r[indx[[1]],indx[[2]]] <- lst[[i]]
106          else ## square matrix          else ## square matrix
107              ret[indx[,1],indx[,2]] <- mlist[[i]]              r[indx[,1], indx[,2]] <- lst[[i]]
108      }      }
109      ## slightly debatable if we really should return Csparse.. :      r
110      as(ret, "CsparseMatrix")  }
111    ## expand(<mer>) needed something like bdiag() for lower-triangular
112    ## (Tsparse) Matrices; hence Doug Bates provided a much more efficient
113    ##  implementation for those; now extended and generalized:
114    .bdiag <- function(lst) {
115        ## block-diagonal matrix [a dgTMatrix] from list of matrices
116        stopifnot(is.list(lst), (nl <- length(lst)) >= 1)
117    
118        Tlst <- lapply(lapply(lst, Matrix:::as_Csp2), # includes "diagU2N"
119                       as, "TsparseMatrix")
120        if(nl == 1) return(Tlst[[1]])
121        ## else
122        i_off <- c(0L, cumsum(vapply(Tlst, nrow, 1L)))
123        j_off <- c(0L, cumsum(vapply(Tlst, ncol, 1L)))
124    
125        clss <- vapply(Tlst, class, "")
126        typ <- substr(clss, 2, 2)
127        knd <- substr(clss, 1, 1)
128        sym <- typ == "s" # symmetric ones
129        tri <- typ == "t" # triangular ones
130        use.n <- any(is.n <- knd == "n")
131        if(use.n && !(use.n <- all(is.n))) {
132            Tlst[is.n] <- lapply(Tlst[is.n], as, "lMatrix")
133            knd [is.n] <- "l"
134        }
135        use.l <- !use.n && all(knd == "l")
136        if(all(sym)) { ## result should be *symmetric*
137            uplos <- vapply(Tlst, slot, ".", "uplo") ## either "U" or "L"
138            tLU <- table(uplos)# of length 1 or 2 ..
139            if(length(tLU) == 1) { ## all "U" or all "L"
140                useU <- uplos[1] == "U"
141            } else { ## length(tLU) == 2, counting "L" and "U"
142                useU <- diff(tLU) >= 0
143                if(useU && (hasL <- tLU[1] > 0))
144                    Tlst[hasL] <- lapply(Tlst[hasL], t)
145                else if(!useU && (hasU <- tLU[2] > 0))
146                    Tlst[hasU] <- lapply(Tlst[hasU], t)
147            }
148            if(use.n) { ## return nsparseMatrix :
149                r <- new("nsTMatrix")
150            } else {
151                r <- new(paste0(if(use.l) "l" else "d", "sTMatrix"))
152                r@x <- unlist(lapply(Tlst, slot, "x"))
153            }
154            r@uplo <- if(useU) "U" else "L"
155        }
156        else if(all(tri) && { ULs <- vapply(Tlst, slot, ".", "uplo")##  "U" or "L"
157                              all(ULs[1L] == ULs[-1L]) } ## all upper or all lower
158           ){ ## *triangular* result
159    
160            if(use.n) { ## return nsparseMatrix :
161                r <- new("ntTMatrix")
162            } else {
163                r <- new(paste0(if(use.l) "l" else "d", "tTMatrix"))
164                r@x <- unlist(lapply(Tlst, slot, "x"))
165            }
166            r@uplo <- ULs[1L]
167        }
168        else {
169            if(any(sym))
170                Tlst[sym] <- lapply(Tlst[sym], as, "generalMatrix")
171            if(use.n) { ## return nsparseMatrix :
172                r <- new("ngTMatrix")
173            } else {
174                r <- new(paste0(if(use.l) "l" else "d", "gTMatrix"))
175                r@x <- unlist(lapply(Tlst, slot, "x"))
176            }
177        }
178        r@Dim <- c(i_off[nl+1], j_off[nl + 1])
179        r@i <- unlist(lapply(1:nl, function(k) Tlst[[k]]@i + i_off[k]))
180        r@j <- unlist(lapply(1:nl, function(k) Tlst[[k]]@j + j_off[k]))
181        r
182    }
183    
184    bdiag <- function(...) {
185        if((nA <- nargs()) == 0) return(new("dgCMatrix"))
186        if(nA == 1 && !is.list(...))
187            return(as(..., "CsparseMatrix"))
188        alis <- if(nA == 1 && is.list(..1)) ..1 else list(...)
189        if(length(alis) == 1)
190            return(as(alis[[1]], "CsparseMatrix"))
191    
192        ## else : two or more arguments
193        as(.bdiag(alis), "CsparseMatrix")
194  }  }
195    
196    
197  .diag2tT <- function(from, uplo = "U", kind = .M.kind(from)) {  .diag2tT <- function(from, uplo = "U", kind = .M.kind(from)) {
198      ## to triangular Tsparse      ## to triangular Tsparse
199      i <- if(from@diag == "U") integer(0) else seq_len(from@Dim[1]) - 1L      i <- if(from@diag == "U") integer(0) else seq_len(from@Dim[1]) - 1L
200      new(paste(kind, "tTMatrix", sep=''),      new(paste0(kind, "tTMatrix"),
201          diag = from@diag, Dim = from@Dim, Dimnames = from@Dimnames,          diag = from@diag, Dim = from@Dim, Dimnames = from@Dimnames,
202          uplo = uplo,          uplo = uplo,
203          x = from@x, # <- ok for diag = "U" and "N" (!)          x = from@x, # <- ok for diag = "U" and "N" (!)
# Line 91  Line 208 
208      ## to symmetric Tsparse      ## to symmetric Tsparse
209      n <- from@Dim[1]      n <- from@Dim[1]
210      i <- seq_len(n) - 1L      i <- seq_len(n) - 1L
211      new(paste(kind, "sTMatrix", sep=''),      new(paste0(kind, "sTMatrix"),
212          Dim = from@Dim, Dimnames = from@Dimnames,          Dim = from@Dim, Dimnames = from@Dimnames,
213          i = i, j = i, uplo = uplo,          i = i, j = i, uplo = uplo,
214          x = if(from@diag == "N") from@x else ## "U"-diag          x = if(from@diag == "N") from@x else ## "U"-diag
# Line 116  Line 233 
233          .diag2tT(d, uplo = if(extends(clx,"triangularMatrix")) x@uplo else "U", kind)          .diag2tT(d, uplo = if(extends(clx,"triangularMatrix")) x@uplo else "U", kind)
234  }  }
235    
236    ## FIXME: should not be needed {when ddi* is dsparse* etc}:
237    setMethod("is.finite", signature(x = "diagonalMatrix"),
238              function(x) is.finite(.diag2tT(x)))
239    setMethod("is.infinite", signature(x = "diagonalMatrix"),
240              function(x) is.infinite(.diag2tT(x)))
241    
242  ## In order to evade method dispatch ambiguity warnings,  ## In order to evade method dispatch ambiguity warnings,
243  ## and because we can save a .M.kind() call, we use this explicit  ## and because we can save a .M.kind() call, we use this explicit
# Line 124  Line 246 
246  ## ddi*:  ## ddi*:
247  diag2tT <- function(from) .diag2tT(from, "U", "d")  diag2tT <- function(from) .diag2tT(from, "U", "d")
248  setAs("ddiMatrix", "triangularMatrix", diag2tT)  setAs("ddiMatrix", "triangularMatrix", diag2tT)
249  setAs("ddiMatrix", "sparseMatrix", diag2tT)  ##_no_longer_ setAs("ddiMatrix", "sparseMatrix", diag2tT)
250  ## needed too (otherwise <dense> -> Tsparse is taken):  ## needed too (otherwise <dense> -> Tsparse is taken):
251  setAs("ddiMatrix", "TsparseMatrix", diag2tT)  setAs("ddiMatrix", "TsparseMatrix", diag2tT)
252    setAs("ddiMatrix", "dsparseMatrix", diag2tT)
253  setAs("ddiMatrix", "CsparseMatrix",  setAs("ddiMatrix", "CsparseMatrix",
254        function(from) as(.diag2tT(from, "U", "d"), "CsparseMatrix"))        function(from) as(.diag2tT(from, "U", "d"), "CsparseMatrix"))
255  setAs("ddiMatrix", "symmetricMatrix",  setAs("ddiMatrix", "symmetricMatrix",
# Line 135  Line 258 
258  ## ldi*:  ## ldi*:
259  diag2tT <- function(from) .diag2tT(from, "U", "l")  diag2tT <- function(from) .diag2tT(from, "U", "l")
260  setAs("ldiMatrix", "triangularMatrix", diag2tT)  setAs("ldiMatrix", "triangularMatrix", diag2tT)
261  setAs("ldiMatrix", "sparseMatrix", diag2tT)  ##_no_longer_ setAs("ldiMatrix", "sparseMatrix", diag2tT)
262  ## needed too (otherwise <dense> -> Tsparse is taken):  ## needed too (otherwise <dense> -> Tsparse is taken):
263  setAs("ldiMatrix", "TsparseMatrix", diag2tT)  setAs("ldiMatrix", "TsparseMatrix", diag2tT)
264    setAs("ldiMatrix", "lsparseMatrix", diag2tT)
265  setAs("ldiMatrix", "CsparseMatrix",  setAs("ldiMatrix", "CsparseMatrix",
266        function(from) as(.diag2tT(from, "U", "l"), "CsparseMatrix"))        function(from) as(.diag2tT(from, "U", "l"), "CsparseMatrix"))
267  setAs("ldiMatrix", "symmetricMatrix",  setAs("ldiMatrix", "symmetricMatrix",
# Line 152  Line 276 
276                Dim = from@Dim, Dimnames = from@Dimnames)                Dim = from@Dim, Dimnames = from@Dimnames)
277        })        })
278    
279    setAs("diagonalMatrix", "nsparseMatrix", function(from) as(from, "nMatrix"))
280    
281  ## Cheap fast substitute for diag() which *does* preserve the mode of x :  ## Cheap fast substitute for diag() which *does* preserve the mode of x :
282  mkDiag <- function(x, n) {  mkDiag <- function(x, n) {
283      y <- matrix(as0(mod=mode(x)), n,n)      y <- matrix(as0(mod=mode(x)), n,n)
284      if (n > 0) y[1L + 0:(n - 1L) * (n + 1L)] <- x      if (n > 0) y[1L + 0:(n - 1L) * (n + 1)] <- x
285      y      y
286  }  }
287    
# Line 173  Line 298 
298                mod.x <- mode(x@x)                mod.x <- mode(x@x)
299                r <- vector(mod.x, length = n^2)                r <- vector(mod.x, length = n^2)
300                if(n)                if(n)
301                    r[1 + 0:(n - 1) * (n + 1)] <-                    r[1 + 0:(n - 1L) * (n + 1)] <-
302                        if(x@diag == "U") as1(mod=mod.x) else x@x                        if(x@diag == "U") as1(mod=mod.x) else x@x
303                r                r
304            })            })
# Line 181  Line 306 
306  setAs("diagonalMatrix", "generalMatrix", # prefer sparse:  setAs("diagonalMatrix", "generalMatrix", # prefer sparse:
307        function(from) as(as(from, "CsparseMatrix"), "generalMatrix"))        function(from) as(as(from, "CsparseMatrix"), "generalMatrix"))
308    
309    setAs("diagonalMatrix", "denseMatrix",
310          function(from) as(as(from, "CsparseMatrix"), "denseMatrix"))
311    
312  .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
313    
314  .diag.2N <- function(m) {  .diag.2N <- function(m) {
# Line 201  Line 329 
329            d <- dim(from)            d <- dim(from)
330            if(d[1] != (n <- d[2])) stop("non-square matrix")            if(d[1] != (n <- d[2])) stop("non-square matrix")
331            if(any(from[row(from) != col(from)] != 0))            if(any(from[row(from) != col(from)] != 0))
332                stop("has non-zero off-diagonal entries")                stop("matrix with non-zero off-diagonals cannot be coerced to diagonalMatrix")
333            x <- diag(from)            x <- diag(from)
334            if(is.logical(x)) {            if(is.logical(x)) {
335                cl <- "ldiMatrix"                cl <- "ldiMatrix"
336                uni <- all(x)                uni <- allTrue(x) ## uni := {is it unit-diagonal ?}
337            } else {            } else {
338                cl <- "ddiMatrix"                cl <- "ddiMatrix"
339                uni <- all(x == 1)                uni <- allTrue(x == 1)
340                storage.mode(x) <- "double"                storage.mode(x) <- "double"
341            } ## TODO: complex            } ## TODO: complex
342            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 353 
353            x <- diag(from)            x <- diag(from)
354            if(is.logical(x)) {            if(is.logical(x)) {
355                cl <- "ldiMatrix"                cl <- "ldiMatrix"
356                uni <- all(x)                uni <- allTrue(x)
357            } else {            } else {
358                cl <- "ddiMatrix"                cl <- "ddiMatrix"
359                uni <- all(x == 1)                uni <- allTrue(x == 1)
360                storage.mode(x) <- "double"                storage.mode(x) <- "double"
361            }            } ## TODO: complex
362            new(cl, Dim = c(n,n), diag = if(uni) "U" else "N",            new(cl, Dim = c(n,n), diag = if(uni) "U" else "N",
363                x = if(uni) x[FALSE] else x)                x = if(uni) x[FALSE] else x)
364        })        })
# Line 240  Line 368 
368            function(x = 1, nrow, ncol) .diag.x(x))            function(x = 1, nrow, ncol) .diag.x(x))
369    
370  subDiag <- function(x, i, j, ..., drop) {  subDiag <- function(x, i, j, ..., drop) {
371      x <- as(x, "sparseMatrix")      x <- as(x, "CsparseMatrix") ## << was "TsparseMatrix" (Csparse is faster now)
372      x <- if(missing(i))      x <- if(missing(i))
373          x[, j, drop=drop]          x[, j, drop=drop]
374      else if(missing(j))      else if(missing(j))
375          x[i, , drop=drop]          if(nargs() == 4) x[i, , drop=drop] else x[i, drop=drop]
376      else      else
377          x[i,j, drop=drop]          x[i,j, drop=drop]
378      if(isS4(x) && isDiagonal(x)) as(x, "diagonalMatrix") else x      if(isS4(x) && isDiagonal(x)) as(x, "diagonalMatrix") else x
# Line 254  Line 382 
382                           j = "index", drop = "logical"), subDiag)                           j = "index", drop = "logical"), subDiag)
383  setMethod("[", signature(x = "diagonalMatrix", i = "index",  setMethod("[", signature(x = "diagonalMatrix", i = "index",
384                          j = "missing", drop = "logical"),                          j = "missing", drop = "logical"),
385            function(x, i, j, ..., drop) subDiag(x, i=i, drop=drop))            function(x, i, j, ..., drop) {
386                  na <- nargs()
387                  Matrix.msg("diag[i,m,l] : nargs()=", na, .M.level = 2)
388                  if(na == 4)
389                       subDiag(x, i=i, , drop=drop)
390                  else subDiag(x, i=i,   drop=drop)
391              })
392  setMethod("[", signature(x = "diagonalMatrix", i = "missing",  setMethod("[", signature(x = "diagonalMatrix", i = "missing",
393                           j = "index", drop = "logical"),                           j = "index", drop = "logical"),
394            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 398 
398  ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch  ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch
399  ## Only(?) current bug:  x[i] <- value  is wrong when  i is *vector*  ## Only(?) current bug:  x[i] <- value  is wrong when  i is *vector*
400  replDiag <- function(x, i, j, ..., value) {  replDiag <- function(x, i, j, ..., value) {
401      x <- as(x, "sparseMatrix")      x <- as(x, "CsparseMatrix")# was "Tsparse.." till 2012-07
402      if(missing(i))      if(missing(i))
403          x[, j] <- value          x[, j] <- value
404      else if(missing(j)) { ##  x[i , ] <- v  *OR*   x[i] <- v      else if(missing(j)) { ##  x[i , ] <- v  *OR*   x[i] <- v
# Line 304  Line 438 
438                                   if(any(value != one | is.na(value))) {                                   if(any(value != one | is.na(value))) {
439                                       x@diag <- "N"                                       x@diag <- "N"
440                                       x@x <- rep.int(one, x@Dim[1])                                       x@x <- rep.int(one, x@Dim[1])
441                                   }                                   } else return(x)
442                               }                               }
443                               x@x[ii] <- value                               x@x[ii] <- value
444                               x                               x
445                           } else { ## no longer diagonal, but remain sparse:                           } else { ## no longer diagonal, but remain sparse:
446                               x <- as(x, "sparseMatrix")                               x <- as(x, "TsparseMatrix")
447                               x[i] <- value                               x[i] <- value
448                               x                               x
449                           }                           }
# Line 325  Line 459 
459                                  j = "index", value = "replValue"),                                  j = "index", value = "replValue"),
460                   function(x,i,j, ..., value) replDiag(x, j=j, value=value))                   function(x,i,j, ..., value) replDiag(x, j=j, value=value))
461    
462    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",
463                                    value = "sparseMatrix"),
464                     function (x, i, j, ..., value)
465                     callGeneric(x=x, , j=j, value = as(value, "sparseVector")))
466    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "missing",
467                                    value = "sparseMatrix"),
468                     function (x, i, j, ..., value)
469                     callGeneric(x=x, i=i, , value = as(value, "sparseVector")))
470    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "index",
471                                    value = "sparseMatrix"),
472                     function (x, i, j, ..., value)
473                     callGeneric(x=x, i=i, j=j, value = as(value, "sparseVector")))
474    
475    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",
476                                    value = "sparseVector"),
477                     replDiag)
478    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "missing",
479                                    value = "sparseVector"),
480                     replDiag)
481    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "index",
482                                    value = "sparseVector"),
483                     replDiag)
484    
485    
486  setMethod("t", signature(x = "diagonalMatrix"),  setMethod("t", signature(x = "diagonalMatrix"),
487            function(x) { x@Dimnames <- x@Dimnames[2:1] ; x })            function(x) { x@Dimnames <- x@Dimnames[2:1] ; x })
# Line 371  Line 528 
528  ##       ---------------------  ##       ---------------------
529  ## Note that "ldi" logical are treated as numeric  ## Note that "ldi" logical are treated as numeric
530  diagdiagprod <- function(x, y) {  diagdiagprod <- function(x, y) {
531      n <- dimCheck(x,y)[1]      dimCheck(x,y)
532      if(x@diag != "U") {      if(x@diag != "U") {
533          if(y@diag != "U") {          if(y@diag != "U") {
534              nx <- x@x * y@x              nx <- x@x * y@x
# Line 403  Line 560 
560      dx <- dim(x)      dx <- dim(x)
561      dy <- dim(y)      dy <- dim(y)
562      if(dx[2] != dy[1]) stop("non-matching dimensions")      if(dx[2] != dy[1]) stop("non-matching dimensions")
     n <- dx[1]  
563      as(if(x@diag == "U") y else x@x * y, "Matrix")      as(if(x@diag == "U") y else x@x * y, "Matrix")
564  }  }
565  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),
# Line 472  Line 628 
628  ##        function(x, y = NULL) {  ##        function(x, y = NULL) {
629  ##           })  ##           })
630    
631  ## setMethod("tcrossprod", signature(x = "denseMatrix", y = "diagonalMatrix"),  Cspdiagprod <- function(x, y) {
632  ##        function(x, y = NULL) {      dx <- dim(x <- .Call(Csparse_diagU2N, x))
633  ##           })      dy <- dim(y)
634        if(dx[2] != dy[1]) stop("non-matching dimensions")
635        if(y@diag == "N") {
636            if(!all(y@x[1L] == y@x[-1L]) && is(x, "symmetricMatrix"))
637                x <- as(x, "generalMatrix")
638            ind <- rep.int(seq_len(dx[2]), x@p[-1] - x@p[-dx[2]-1L])
639            x@x <- x@x * y@x[ind]
640        }
641        if(is(x, "compMatrix") && length(xf <- x@factors)) {
642            ## instead of dropping all factors, be smart about some
643            ## TODO ......
644            x@factors <- list()
645        }
646        x
647    }
648    
649    diagCspprod <- function(x, y) {
650        dx <- dim(x)
651        dy <- dim(y <- .Call(Csparse_diagU2N, y))
652        if(dx[2] != dy[1]) stop("non-matching dimensions")
653        if(x@diag == "N") {
654            if(!all(x@x[1L] == x@x[-1L]) && is(y, "symmetricMatrix"))
655                y <- as(y, "generalMatrix")
656            y@x <- y@x * x@x[y@i + 1L]
657        }
658        if(is(y, "compMatrix") && length(yf <- y@factors)) {
659            ## instead of dropping all factors, be smart about some
660            ## TODO
661            keep <- character()
662            if(iLU <- names(yf) == "LU") {
663                ## TODO keep <- "LU"
664            }
665            y@factors <- yf[keep]
666        }
667        y
668    }
669    
670    setMethod("crossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
671              function(x, y = NULL) diagCspprod(x, y))
672    
673  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
674            function(x, y = NULL) crossprod(as(x, "sparseMatrix"), y))            function(x, y = NULL) diagCspprod(x, as(y, "CsparseMatrix")))
675    
676    ## Prefer calling diagCspprod to Cspdiagprod if going to transpose anyway
677    ##  x'y == (y'x)'
678    setMethod("crossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
679              function(x, y = NULL) t(diagCspprod(y, x)))
680    
681  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
682            function(x, y = NULL) crossprod(x, as(y, "sparseMatrix")))            function(x, y = NULL) t(diagCspprod(y, as(x, "Csparsematrix"))))
683    
684    setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
685              function(x, y = NULL) diagCspprod(x, t(y)))
686    
687  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
688            function(x, y = NULL) tcrossprod(as(x, "sparseMatrix"), y))            function(x, y = NULL) diagCspprod(x, t(as(y, "CsparseMatrix"))))
689    
690    setMethod("tcrossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
691              function(x, y = NULL) Cspdiagprod(x, y))
692    
693  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
694            function(x, y = NULL) tcrossprod(x, as(y, "sparseMatrix")))            function(x, y = NULL) Cspdiagprod(as(x, "CsparseMatrix"), y))
695    
696    setMethod("%*%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
697              function(x, y) diagCspprod(x, y))
698    
 ## FIXME?: In theory, this can be done *FASTER*, in some cases, via tapply1()  
699  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
700            function(x, y) as(x, "sparseMatrix") %*% y)            function(x, y) diagCspprod(as(x, "CsparseMatrix"), y))
701    
702  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
703            function(x, y) x %*% as(y, "sparseMatrix"))            function(x, y) Cspdiagprod(as(x, "CsparseMatrix"), y))
704  ## 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)  
705  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
706            function(x, y) x %*% as(y, "CsparseMatrix"))            function(x, y) Cspdiagprod(x, y))
707  ## NB: this is *not* needed for Tsparse & Rsparse  
708  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*
709  ##       do indeed work by going through sparse (and *not* ddense)!  ##       do indeed work by going through sparse (and *not* ddense)!
710    
   
   
711  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),
712            function(a, b, ...) {            function(a, b, ...) {
713                a@x <- 1/ a@x                a@x <- 1/ a@x
# Line 515  Line 716 
716            })            })
717    
718  solveDiag <- function(a, b, ...) {  solveDiag <- function(a, b, ...) {
719      if((n <- a@Dim[1]) != nrow(b))      if(a@Dim[1] != nrow(b))
720          stop("incompatible matrix dimensions")          stop("incompatible matrix dimensions")
721      ## trivially invert a 'in place' and multiply:      ## trivially invert a 'in place' and multiply:
722      a@x <- 1/ a@x      a@x <- 1/ a@x
# Line 534  Line 735 
735  ###---------------- <Ops> (<Arith>, <Logic>, <Compare> ) ----------------------  ###---------------- <Ops> (<Arith>, <Logic>, <Compare> ) ----------------------
736    
737  ## Use function for several signatures, in order to evade  ## Use function for several signatures, in order to evade
 ## ambiguous dispatch for "ddi", since there's also Arith(ddense., ddense.)  
738  diagOdiag <- function(e1,e2) {  diagOdiag <- function(e1,e2) {
739      ## result should also be diagonal _ if possible _      ## result should also be diagonal _ if possible _
740      r <- callGeneric(.diag.x(e1), .diag.x(e2)) # error if not "compatible"      r <- callGeneric(.diag.x(e1), .diag.x(e2)) # error if not "compatible"
# Line 542  Line 742 
742      r00 <- callGeneric(if(is.numeric(e1@x)) 0 else FALSE,      r00 <- callGeneric(if(is.numeric(e1@x)) 0 else FALSE,
743                         if(is.numeric(e2@x)) 0 else FALSE)                         if(is.numeric(e2@x)) 0 else FALSE)
744      if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* diagonal      if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* diagonal
745          if(is.numeric(r)) {          if(is.numeric(r)) { # "double" *or* "integer"
746              if(is.numeric(e2@x)) {              if(is.numeric(e2@x)) {
747                  e2@x <- r; return(.diag.2N(e2)) }                  e2@x <- r; return(.diag.2N(e2)) }
748              if(!is.numeric(e1@x))              if(!is.numeric(e1@x))
749                  ## e.g. e1, e2 are logical;                  ## e.g. e1, e2 are logical;
750                  e1 <- as(e1, "dMatrix")                  e1 <- as(e1, "dMatrix")
751                if(!is.double(r)) r <- as.double(r)
752          }          }
753          else if(is.logical(r))          else if(is.logical(r))
754              e1 <- as(e1, "lMatrix")              e1 <- as(e1, "lMatrix")
# Line 556  Line 757 
757          .diag.2N(e1)          .diag.2N(e1)
758      }      }
759      else { ## result not diagonal, but at least symmetric:      else { ## result not diagonal, but at least symmetric:
760            ## e.g., m == m
761          isNum <- (is.numeric(r) || is.numeric(r00))          isNum <- (is.numeric(r) || is.numeric(r00))
762          isLog <- (is.logical(r) || is.logical(r00))          isLog <- (is.logical(r) || is.logical(r00))
763            Matrix.msg("exploding <diag> o <diag> into dense matrix", .M.level = 2)
         if(getOption("verbose"))  
             message("exploding  <diag>  o  <diag>  into dense matrix")  
764          d <- e1@Dim          d <- e1@Dim
765          n <- d[1]          n <- d[1]
766          stopifnot(length(r) == n)          stopifnot(length(r) == n)
767            if(isNum && !is.double(r)) r <- as.double(r)
768          xx <- as.vector(matrix(rbind(r, matrix(r00,n,n)), n,n))          xx <- as.vector(matrix(rbind(r, matrix(r00,n,n)), n,n))
769          newcl <-          newcl <-
770              paste(if(isNum) "d" else if(isLog) {              paste0(if(isNum) "d" else if(isLog) {
771                  if(!any(is.na(r)) && !any(is.na(r00))) "n" else "l"                  if(!any(is.na(r)) && !any(is.na(r00))) "n" else "l"
772              } else stop("not yet implemented .. please report")              } else stop("not yet implemented .. please report"), "syMatrix")
                   ,  
                   "syMatrix", sep='')  
773    
774          new(newcl, Dim = e1@Dim, Dimnames = e1@Dimnames, x = xx)          new(newcl, Dim = e1@Dim, Dimnames = e1@Dimnames, x = xx)
775      }      }
# Line 588  Line 787 
787              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)
788  }  }
789    
790  ## FIXME:    diagonal  o  triangular  |-->  triangular  ## diagonal  o  triangular  |-->  triangular
791  ## -----     diagonal  o  symmetric   |-->  symmetric  ## diagonal  o  symmetric   |-->  symmetric
792  ##    {also when other is sparse: do these "here" --  ##    {also when other is sparse: do these "here" --
793  ##     before conversion to sparse, since that loses "diagonality"}  ##     before conversion to sparse, since that loses "diagonality"}
794    diagOtri <- function(e1,e2) {
795        ## result must be triangular
796        r <- callGeneric(d1 <- .diag.x(e1), diag(e2)) # error if not "compatible"
797        ## Check what happens with non-diagonals, i.e. (0 o 0), (FALSE o 0), ...:
798        e1.0 <- if(.n1 <- is.numeric(d1   )) 0 else FALSE
799        r00 <- callGeneric(e1.0, if(.n2 <- is.numeric(e2[0])) 0 else FALSE)
800        if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* triangular
801            diag(e2) <- r
802            ## check what happens "in the triangle"
803            e2.2 <- if(.n2) 2 else TRUE
804            if(!callGeneric(e1.0, e2.2) == e2.2) { # values "in triangle" can change:
805                n <- dim(e2)[1L]
806                it <- indTri(n, upper = (e2@uplo == "U"))
807                e2[it] <- callGeneric(e1.0, e2[it])
808            }
809            e2
810        }
811        else { ## result not triangular ---> general
812            rr <- as(e2, "generalMatrix")
813            diag(rr) <- r
814            rr
815        }
816    }
817    
818    
819    setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "triangularMatrix"),
820              diagOtri)
821    ## For the reverse,  Ops == "Arith" | "Compare" | "Logic"
822    ##   'Arith'  :=  '"+"', '"-"', '"*"', '"^"', '"%%"', '"%/%"', '"/"'
823    setMethod("Arith", signature(e1 = "triangularMatrix", e2 = "diagonalMatrix"),
824              function(e1,e2)
825          { ## this must only trigger for *dense* e1
826              switch(.Generic,
827                     "+" = .Call(dtrMatrix_addDiag, as(e1,"dtrMatrix"),   .diag.x(e2)),
828                     "-" = .Call(dtrMatrix_addDiag, as(e1,"dtrMatrix"), - .diag.x(e2)),
829                     "*" = {
830                         n <- e2@Dim[1L]
831                         d2 <- if(e2@diag == "U") { # unit-diagonal
832                             d <- rep.int(as1(e2@x), n)
833                             e2@x <- d
834                             e2@diag <- "N"
835                             d
836                         } else e2@x
837                         e2@x <- diag(e1) * d2
838                         e2
839                     },
840                     "^" = { ## will be dense ( as  <ANY> ^ 0 == 1 ):
841                         e1 ^ as(e2, "denseMatrix")
842                     },
843                     ## otherwise:
844                     callGeneric(e1, diag2Tsmart(e2,e1)))
845    })
846    
847    ## Compare --> 'swap' (e.g.   e1 < e2   <==>  e2 > e1 ):
848    setMethod("Compare", signature(e1 = "triangularMatrix", e2 = "diagonalMatrix"),
849              .Cmp.swap)
850    ## '&' and "|'  are commutative:
851    setMethod("Logic", signature(e1 = "triangularMatrix", e2 = "diagonalMatrix"),
852              function(e1,e2) callGeneric(e2, e1))
853    
854  ## For almost everything else, diag* shall be treated "as sparse" :  ## For almost everything else, diag* shall be treated "as sparse" :
855  ## These are cheap implementations via coercion  ## These are cheap implementations via coercion
# Line 615  Line 873 
873  ##       Compare --> logical  ##       Compare --> logical
874  ##       Logic   --> logical: "lMatrix"  ##       Logic   --> logical: "lMatrix"
875    
876  ##  other = "numeric" : stay diagonal if possible  ## Other = "numeric" : stay diagonal if possible
877  ## ddi*: Arith: result numeric, potentially ddiMatrix  ## ddi*: Arith: result numeric, potentially ddiMatrix
878  setMethod("Arith", signature(e1 = "ddiMatrix", e2 = "numeric"),  for(arg2 in c("numeric","logical"))
879    setMethod("Arith", signature(e1 = "ddiMatrix", e2 = arg2),
880            function(e1,e2) {            function(e1,e2) {
881                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]
882                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
883                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
884                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
885                    if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)                    if(e1@diag == "U") {
886                    if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {                        if(any((r <- callGeneric(1, e2)) != 1)) {
887                        e1@diag <- "N"                        e1@diag <- "N"
888                        if(L1) r <- rep.int(r, n)                            e1@x[seq_len(n)] <- r # possibly recycling r
889                    } else                        } ## else: result = e1  (is "U" diag)
890                      } else {
891                        r <- callGeneric(e1@x, e2)                        r <- callGeneric(e1@x, e2)
892                    e1@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]                        ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
893                    return(e1)                        e1@x[] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
894                }                }
895                      e1
896                  } else
897                callGeneric(diag2tT.u(e1,e2, "d"), e2)                callGeneric(diag2tT.u(e1,e2, "d"), e2)
898            })            })
899    
900  setMethod("Arith", signature(e1 = "numeric", e2 = "ddiMatrix"),  for(arg1 in c("numeric","logical"))
901    setMethod("Arith", signature(e1 = arg1, e2 = "ddiMatrix"),
902            function(e1,e2) {            function(e1,e2) {
903                n <- e2@Dim[1]; nsq <- n*n                n <- e2@Dim[1]
904                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
905                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
906                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
907                    if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)                    if(e2@diag == "U") {
908                    if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {                        if(any((r <- callGeneric(e1, 1)) != 1)) {
909                        e2@diag <- "N"                        e2@diag <- "N"
910                        if(L1) r <- rep.int(r, n)                            e2@x[seq_len(n)] <- r # possibly recycling r
911                    } else                        } ## else: result = e2  (is "U" diag)
912                      } else {
913                        r <- callGeneric(e1, e2@x)                        r <- callGeneric(e1, e2@x)
914                    e2@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]                        ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
915                    return(e2)                        e2@x[] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
916                }                }
917                      e2
918                  } else
919                callGeneric(e1, diag2tT.u(e2,e1, "d"))                callGeneric(e1, diag2tT.u(e2,e1, "d"))
920            })            })
921    
922  ## ldi* Arith --> result numeric, potentially ddiMatrix  ## ldi* Arith --> result numeric, potentially ddiMatrix
923  setMethod("Arith", signature(e1 = "ldiMatrix", e2 = "numeric"),  for(arg2 in c("numeric","logical"))
924    setMethod("Arith", signature(e1 = "ldiMatrix", e2 = arg2),
925            function(e1,e2) {            function(e1,e2) {
926                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]
927                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
928                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
929                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
930                    if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)                    E <- copyClass(e1, "ddiMatrix", c("diag", "Dim", "Dimnames"))#FIXME: if ok, check=FALSE
931                    if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {                    ## E <- copyClass(e1, "ddiMatrix", check=FALSE)
932                        e1@diag <- "N"                    ## storage.mode(E@x) <- "double"
933                        if(L1) r <- rep.int(r, n)                    if(e1@diag == "U") {
934                    } else                        if(any((r <- callGeneric(1, e2)) != 1)) {
935                              E@diag <- "N"
936                              E@x[seq_len(n)] <- r # possibly recycling r
937                          } ## else: result = E  (is "U" diag)
938                      } else {
939                        r <- callGeneric(e1@x, e2)                        r <- callGeneric(e1@x, e2)
940                    e1 <- copyClass(e1, "ddiMatrix", c("diag", "Dim", "Dimnames"))                        ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
941                    e1@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]                        E@x[seq_len(n)] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
                   return(e1)  
942                }                }
943                callGeneric(diag2tT.u(e1,e2, "d"), e2)                    E
944                  } else
945                      callGeneric(diag2tT.u(e1,e2, "l"), e2)
946            })            })
947    
948  setMethod("Arith", signature(e1 = "numeric", e2 = "ldiMatrix"),  for(arg1 in c("numeric","logical"))
949    setMethod("Arith", signature(e1 = arg1, e2 = "ldiMatrix"),
950            function(e1,e2) {            function(e1,e2) {
951                n <- e2@Dim[1]; nsq <- n*n                n <- e2@Dim[1]
952                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
953                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
954                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
955                    if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)                    E <- copyClass(e2, "ddiMatrix", c("diag", "Dim", "Dimnames"))#FIXME: if ok, check=FALSE
956                    if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {                    ## E <- copyClass(e2, "ddiMatrix", check=FALSE)
957                        e2@diag <- "N"                    ## storage.mode(E@x) <- "double"
958                        if(L1) r <- rep.int(r, n)                    if(e2@diag == "U") {
959                    } else                        if(any((r <- callGeneric(e1, 1)) != 1)) {
960                              E@diag <- "N"
961                              E@x[seq_len(n)] <- r # possibly recycling r
962                          } ## else: result = E  (is "U" diag)
963                      } else {
964                        r <- callGeneric(e1, e2@x)                        r <- callGeneric(e1, e2@x)
965                    e2 <- copyClass(e2, "ddiMatrix", c("diag", "Dim", "Dimnames"))                        ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
966                    e2@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]                        E@x[seq_len(n)] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
                   return(e2)  
967                }                }
968                callGeneric(e1, diag2tT.u(e2,e1, "d"))                    E
969                  } else
970                      callGeneric(e1, diag2tT.u(e2,e1, "l"))
971            })            })
972    
973  ## ddi*: for "Ops" without Arith --> result logical, potentially ldi  ## ddi*: for "Ops" without "Arith": <Compare> or <Logic> --> result logical, potentially ldi
974  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "numeric"),  ##
975    ## Note that  ("numeric", "ddiMatrix")  is simply swapped, e.g.,
976    if(FALSE) {
977        selectMethod("<", c("numeric","lMatrix"))# Compare
978        selectMethod("&", c("numeric","lMatrix"))# Logic
979    } ## so we don't need to define a method here :
980    
981    for(arg2 in c("numeric","logical"))
982    setMethod("Ops", signature(e1 = "ddiMatrix", e2 = arg2),
983            function(e1,e2) {            function(e1,e2) {
984                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]
985                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
986                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
987                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
988                    if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)                    E <- copyClass(e1, "ldiMatrix", c("diag", "Dim", "Dimnames"))#FIXME: if ok, check=FALSE
989                    if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {                    ## E <- copyClass(e1, "ldiMatrix", check=FALSE)
990                        e1@diag <- "N"                    ## storage.mode(E@x) <- "logical"
991                        if(L1) r <- rep.int(r, n)                    if(e1@diag == "U") {
992                    } else                        if(any((r <- callGeneric(1, e2)) != 1)) {
993                              E@diag <- "N"
994                              E@x[seq_len(n)] <- r # possibly recycling r
995                          } ## else: result = E  (is "U" diag)
996                      } else {
997                        r <- callGeneric(e1@x, e2)                        r <- callGeneric(e1@x, e2)
998                    e1 <- copyClass(e1, "ldiMatrix", c("diag", "Dim", "Dimnames"))                        ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
999                    e1@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]                        E@x[seq_len(n)] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
                   return(e1)  
1000                }                }
1001                callGeneric(diag2tT.u(e1,e2, "l"), e2)                    E
           })  
   
 setMethod("Ops", signature(e1 = "numeric", e2 = "ddiMatrix"),  
           function(e1,e2) {  
               n <- e2@Dim[1]; nsq <- n*n  
               f0 <- callGeneric(e1, 0)  
               if(all(is0(f0))) { # remain diagonal  
                   L1 <- (le <- length(e1)) == 1L  
                   if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)  
                   if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {  
                       e2@diag <- "N"  
                       if(L1) r <- rep.int(r, n)  
1002                    } else                    } else
1003                        r <- callGeneric(e1, e2@x)                    callGeneric(diag2tT.u(e1,e2, "d"), e2)
                   e2 <- copyClass(e2, "ldiMatrix", c("diag", "Dim", "Dimnames"))  
                   e2@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]  
                   return(e2)  
               }  
               callGeneric(e1, diag2tT.u(e2,e1, "l"))  
1004            })            })
1005    
1006  ## ldi*: for "Ops" without Arith --> result logical, potentially ldi  ## ldi*: for "Ops" without "Arith": <Compare> or <Logic> --> result logical, potentially ldi
1007  setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "numeric"),  for(arg2 in c("numeric","logical"))
1008    setMethod("Ops", signature(e1 = "ldiMatrix", e2 = arg2),
1009            function(e1,e2) {            function(e1,e2) {
1010                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]
1011                f0 <- callGeneric(FALSE, e2)                f0 <- callGeneric(FALSE, e2)
1012                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
1013                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
1014                    if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)  
1015                    if(e1@diag == "U" && any((r <- callGeneric(TRUE, e2)) != 1)) {                    if(e1@diag == "U") {
1016                          if(any((r <- callGeneric(TRUE, e2)) != 1)) {
1017                        e1@diag <- "N"                        e1@diag <- "N"
1018                        if(L1) r <- rep.int(r, n)                            e1@x[seq_len(n)] <- r # possibly recycling r
1019                    } else                        } ## else: result = e1  (is "U" diag)
1020                      } else {
1021                        r <- callGeneric(e1@x, e2)                        r <- callGeneric(e1@x, e2)
1022                    e1@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]                        ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
1023                    return(e1)                        e1@x[] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
1024                }                }
1025                callGeneric(diag2tT.u(e1,e2, "l"), e2)                    e1
           })  
   
 setMethod("Ops", signature(e1 = "numeric", e2 = "ldiMatrix"),  
           function(e1,e2) {  
               n <- e2@Dim[1]; nsq <- n*n  
               f0 <- callGeneric(e1, FALSE)  
               if(all(is0(f0))) { # remain diagonal  
                   L1 <- (le <- length(e1)) == 1L  
                   if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)  
                   if(e2@diag == "U" && any((r <- callGeneric(e1, TRUE)) != 1)) {  
                       e2@diag <- "N"  
                       if(L1) r <- rep.int(r, n)  
1026                    } else                    } else
1027                        r <- callGeneric(e1, e2@x)                    callGeneric(diag2tT.u(e1,e2, "l"), e2)
                   e2@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]  
                   return(e2)  
               }  
               callGeneric(e1, diag2tT.u(e2,e1, "l"))  
1028            })            })
1029    
1030    
   
1031  ## Not {"sparseMatrix", "numeric} :  {"denseMatrix", "matrix", ... }  ## Not {"sparseMatrix", "numeric} :  {"denseMatrix", "matrix", ... }
1032  for(other in c("ANY", "Matrix", "dMatrix")) {  for(other in c("ANY", "Matrix", "dMatrix")) {
1033      ## ddi*:      ## ddi*:
1034      setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other),      setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other),
1035                function(e1,e2) callGeneric(diag2tT.u(e1,e2, "d"), e2))                function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2))
1036      setMethod("Ops", signature(e1 = other, e2 = "ddiMatrix"),      setMethod("Ops", signature(e1 = other, e2 = "ddiMatrix"),
1037                function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "d")))                function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "d")))
1038      ## ldi*:      ## ldi*:
1039      setMethod("Ops", signature(e1 = "ldiMatrix", e2 = other),      setMethod("Ops", signature(e1 = "ldiMatrix", e2 = other),
1040                function(e1,e2) callGeneric(diag2tT.u(e1,e2, "l"), e2))                function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2))
1041      setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),      setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),
1042                function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "l")))                function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "l")))
1043  }  }
1044  ## This should *not* dispatch to <dense> methods (in ./Ops.R ), as  
1045  ##  FALSE & <anything> |-> FALSE : hence result should be diagonal:  ## Direct subclasses of "denseMatrix": currently ddenseMatrix, ldense... :
1046  for(cl in diCls) {  dense.subCl <- local({ dM.scl <- getClass("denseMatrix")@subclasses
1047      setMethod("&", signature(e1 = cl, e2 = "ANY"),                         names(dM.scl)[vapply(dM.scl, slot, 0, "distance") == 1] })
1048                function(e1,e2) e1 & as(e2,"Matrix"))  for(DI in diCls) {
1049      setMethod("&", signature(e1 = "ANY", e2 = cl),      dMeth <- if(extends(DI, "dMatrix"))
1050                function(e1,e2) as(e1,"Matrix") & e2)          function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2)
1051      for(c2 in c("denseMatrix", "Matrix")) {      else # "lMatrix", the only other kind for now
1052          setMethod("&", signature(e1 = cl, e2 = c2),          function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2)
1053                    function(e1,e2) e1 & Diagonal(x = diag(e2)))      for(c2 in c(dense.subCl, "Matrix")) {
1054          setMethod("&", signature(e1 = c2, e2 = cl),          for(Fun in c("*", "&")) {
1055                    function(e1,e2) Diagonal(x = diag(e1)) & e2)              setMethod(Fun, signature(e1 = DI, e2 = c2),
1056                          function(e1,e2) callGeneric(e1, Diagonal(x = diag(e2))))
1057                setMethod(Fun, signature(e1 = c2, e2 = DI),
1058                          function(e1,e2) callGeneric(Diagonal(x = diag(e1)), e2))
1059            }
1060            setMethod("^", signature(e1 = c2, e2 = DI),
1061                      function(e1,e2) callGeneric(Diagonal(x = diag(e1)), e2))
1062            for(Fun in c("%%", "%/%", "/")) ## 0 <op> 0 |--> NaN  for these.
1063                setMethod(Fun, signature(e1 = DI, e2 = c2), dMeth)
1064      }      }
1065  }  }
1066    
1067    
1068  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"
1069  ### ----------  any, all: separately here  ### ----------   the last 4: separately here
1070  for(cl in diCls) {  for(cl in diCls) {
1071  setMethod("any", cl,  setMethod("any", cl,
1072            function (x, ..., na.rm) {            function (x, ..., na.rm) {
1073                if(any(x@Dim == 0)) FALSE                if(any(x@Dim == 0)) FALSE
1074                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)
1075            })            })
1076  setMethod("all",  cl, function (x, ..., na.rm) any(x@Dim == 0))  setMethod("all",  cl, function (x, ..., na.rm) {
1077  setMethod("prod", cl, function (x, ..., na.rm) as.numeric(any(x@Dim == 0)))      n <- x@Dim[1]
1078        if(n >= 2) FALSE
1079        else if(n == 0 || x@diag == "U") TRUE
1080        else all(x@x, ..., na.rm = na.rm)
1081    })
1082    setMethod("prod", cl, function (x, ..., na.rm) {
1083        n <- x@Dim[1]
1084        if(n >= 2) 0
1085        else if(n == 0 || x@diag == "U") 1
1086        else ## n == 1, diag = "N" :
1087            prod(x@x, ..., na.rm = na.rm)
1088    })
1089    
1090  setMethod("sum", cl,  setMethod("sum", cl,
1091            function(x, ..., na.rm) {            function(x, ..., na.rm) {
# Line 847  Line 1124 
1124      invisible(x)      invisible(x)
1125  }  }
1126    
1127    ## somewhat consistent with "print" for sparseMatrix :
1128    setMethod("print", signature(x = "diagonalMatrix"), prDiag)
1129    
1130  setMethod("show", signature(object = "diagonalMatrix"),  setMethod("show", signature(object = "diagonalMatrix"),
1131            function(object) {            function(object) {
1132                d <- dim(object)                d <- dim(object)
1133                cl <- class(object)                cl <- class(object)
1134                cat(sprintf('%d x %d diagonal matrix of class "%s"\n',                cat(sprintf('%d x %d diagonal matrix of class "%s"',
1135                            d[1], d[2], cl))                            d[1], d[2], cl))
1136                  if(d[1] < 50) {
1137                      cat("\n")
1138                prDiag(object)                prDiag(object)
1139                  } else {
1140                      cat(", with diagonal entries\n")
1141                      show(diag(object))
1142                      invisible(object)
1143                  }
1144            })            })
1145    
1146    rm(dense.subCl, diCls)# not used elsewhere
1147    
1148    setMethod("summary", signature(object = "diagonalMatrix"),
1149              function(object, ...) {
1150                  d <- dim(object)
1151                  r <- summary(object@x, ...)
1152                  attr(r, "header") <-
1153                      sprintf('%d x %d diagonal Matrix of class "%s"',
1154                              d[1], d[2], class(object))
1155                  ## use ole' S3 technology for such a simple case
1156                  class(r) <- c("diagSummary", class(r))
1157                  r
1158              })
1159    
1160    print.diagSummary <- function (x, ...) {
1161        cat(attr(x, "header"),"\n")
1162        class(x) <- class(x)[-1]
1163        print(x, ...)
1164        invisible(x)
1165    }

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

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