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 2237, Fri Jul 25 06:55:42 2008 UTC pkg/Matrix/R/diagMatrix.R revision 2984, Sat Apr 12 21:37:37 2014 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))      n <- if(missing(n)) length(x) else {
         n <- length(x)  
     else {  
10          stopifnot(length(n) == 1, n == as.integer(n), n >= 0)          stopifnot(length(n) == 1, n == as.integer(n), n >= 0)
11          n <- as.integer(n)          as.integer(n)
12      }      }
13    
14      if(missing(x)) ## unit diagonal matrix      if(missing(x)) ## unit diagonal matrix
15          new("ddiMatrix", Dim = c(n,n), diag = "U")          new("ddiMatrix", Dim = c(n,n), diag = "U")
16      else {      else {
17          lx <- length(x)          lx <- length(x)
18          stopifnot(lx == 1 || lx == n) # but keep 'x' short for now          lx.1 <- lx == 1L
19            stopifnot(lx.1 || lx == n) # but keep 'x' short for now
20          if(is.logical(x))          if(is.logical(x))
21              cl <- "ldiMatrix"              cl <- "ldiMatrix"
22          else if(is.numeric(x)) {          else if(is.numeric(x)) {
# Line 27  Line 26 
26          else if(is.complex(x)) {          else if(is.complex(x)) {
27              cl <- "zdiMatrix"  # will not yet work              cl <- "zdiMatrix"  # will not yet work
28          } else stop("'x' has invalid data type")          } else stop("'x' has invalid data type")
29            if(lx.1 && !is.na(x) && x == 1) # cheap check for uni-diagonal..
30                new(cl, Dim = c(n,n), diag = "U")
31            else
32          new(cl, Dim = c(n,n), diag = "N",          new(cl, Dim = c(n,n), diag = "N",
33              x = if(lx == 1) rep.int(x,n) else x)                  x = if(lx.1) rep.int(x,n) else x)
34      }      }
35  }  }
36    
37  .sparseDiagonal <- function(n, x = rep.int(1,n), uplo = "U", shape = "t") {  .sparseDiagonal <- function(n, x = 1, uplo = "U",
38                                shape = if(missing(cols)) "t" else "g",
39                                unitri, kind,
40                                cols = if(n) 0:(n - 1L) else integer(0))
41    {
42      stopifnot(n == (n. <- as.integer(n)), (n <- n.) >= 0)      stopifnot(n == (n. <- as.integer(n)), (n <- n.) >= 0)
43      if((lx <- length(x)) == 1) x <- rep.int(x, n)      if(!(mcols <- missing(cols)))
44      else if(lx != n) stop("length(x) must be 1 or n")          stopifnot(0 <= (cols <- as.integer(cols)), cols < n)
45      stopifnot(is.character(shape), nchar(shape) == 1,      m <- length(cols)
46                any(shape == c("t","s","g"))) # triangular / symmetric / general      if(missing(kind))
47      kind <-      kind <-
48          if(is.double(x)) "d"          if(is.double(x)) "d"
49          else if(is.logical(x)) "l"          else if(is.logical(x)) "l"
# Line 45  Line 51 
51              storage.mode(x) <- "double"              storage.mode(x) <- "double"
52              "d"              "d"
53          }          }
54      new(paste(kind, shape, "CMatrix", sep=''),      else stopifnot(any(kind == c("d","l","n")))
55          Dim = c(n,n), x = x, uplo = uplo,      stopifnot(is.character(shape), nchar(shape) == 1,
56          i = if(n) 0:(n - 1L) else integer(0), p = 0:n)                any(shape == c("t","s","g"))) # triangular / symmetric / general
57        if((missing(unitri) || unitri) && shape == "t" &&
58           (mcols || cols == 0:(n-1L)) &&
59           ((any(kind == c("l", "n")) && allTrue(x)) ||
60            (    kind == "d"          && allTrue(x == 1)))) { ## uni-triangular
61            new(paste0(kind,"tCMatrix"), Dim = c(n,n),
62                       uplo = uplo, diag = "U", p = rep.int(0L, n+1L))
63        }
64        else if(kind == "n") {
65            if(shape == "g")
66                new("ngCMatrix", Dim = c(n,m), i = cols, p = 0:m)
67            else new(paste0("n", shape, "CMatrix"), Dim = c(n,m), uplo = uplo,
68                     i = cols, p = 0:m)
69        }
70        else { ## kind != "n" -- have x slot :
71            if((lx <- length(x)) == 1) x <- rep.int(x, m)
72            else if(lx != m) stop("length(x) must be either 1 or #{cols}")
73            if(shape == "g")
74                new(paste0(kind, "gCMatrix"), Dim = c(n,m),
75                    x = x, i = cols, p = 0:m)
76            else new(paste0(kind, shape, "CMatrix"), Dim = c(n,m), uplo = uplo,
77                     x = x, i = cols, p = 0:m)
78        }
79  }  }
80    
81  ## Pkg 'spdep' had (relatively slow) versions of this as_dsCMatrix_I()  ## Pkg 'spdep' had (relatively slow) versions of this as_dsCMatrix_I()
82  .symDiagonal <- function(n, x = rep.int(1,n), uplo = "U")  .symDiagonal <- function(n, x = rep.int(1,n), uplo = "U")
83      .sparseDiagonal(n, x, uplo, shape = "s")      .sparseDiagonal(n, x, uplo, shape = "s")
84    
85  ## instead of   diagU2N(as(Diagonal(n), "CsparseMatrix")), diag = "N" in any case:  # instead of   diagU2N(as(Diagonal(n), "CsparseMatrix")), diag = "N" in any case:
86  .trDiagonal <- function(n, x = rep.int(1,n), uplo = "U")  .trDiagonal <- function(n, x = 1, uplo = "U", unitri=TRUE)
87      .sparseDiagonal(n, x, uplo, shape = "t")      .sparseDiagonal(n, x, uplo, shape = "t", unitri=unitri)
88    
89    
90  ### 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.
91  ### 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
92  ### 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
93  ### who posted his bdiag() function written in December 1995.  ## who posted his bdiag() function written in December 1995.
94  if(FALSE)##--- no longer used:  if(FALSE)##--- no longer used:
95  .bdiag <- function(lst) {  .bdiag <- function(lst) {
96      ### block-diagonal matrix [a dgTMatrix] from list of matrices      ## block-diagonal matrix [a dgTMatrix] from list of matrices
97      stopifnot(is.list(lst), length(lst) >= 1)      stopifnot(is.list(lst), length(lst) >= 1)
98      dims <- sapply(lst, dim, USE.NAMES=FALSE)      dims <- vapply(lst, dim, 1L, USE.NAMES=FALSE)
99      ## make sure we had all matrices:      ## make sure we had all matrices:
100      if(!(is.matrix(dims) && nrow(dims) == 2))      if(!(is.matrix(dims) && nrow(dims) == 2))
101          stop("some arguments are not matrices")          stop("some arguments are not matrices")
# Line 92  Line 120 
120      ## block-diagonal matrix [a dgTMatrix] from list of matrices      ## block-diagonal matrix [a dgTMatrix] from list of matrices
121      stopifnot(is.list(lst), (nl <- length(lst)) >= 1)      stopifnot(is.list(lst), (nl <- length(lst)) >= 1)
122    
123      Tlst <- lapply(lapply(lst, Matrix:::as_Csp2), # includes "diagU2N"      Tlst <- lapply(lapply(lst, as_Csp2), # includes "diagU2N"
124                     as, "TsparseMatrix")                     as, "TsparseMatrix")
125      if(nl == 1) return(Tlst[[1]])      if(nl == 1) return(Tlst[[1]])
126      ## else      ## else
127      i_off <- c(0L, cumsum(sapply(Tlst, nrow)))      i_off <- c(0L, cumsum(vapply(Tlst, nrow, 1L)))
128      j_off <- c(0L, cumsum(sapply(Tlst, ncol)))      j_off <- c(0L, cumsum(vapply(Tlst, ncol, 1L)))
129    
130      clss <- sapply(Tlst, class)      clss <- vapply(Tlst, class, "")
131      knds <- substr(clss, 2, 2)      typ <- substr(clss, 2, 2)
132      sym  <- knds == "s" # symmetric ones      knd <- substr(clss, 1, 1)
133      tri  <- knds == "t" # triangular ones      sym <- typ == "s" # symmetric ones
134      use.n <- any(is.n <- substr(clss,1,1) == "n")      tri <- typ == "t" # triangular ones
135      if(use.n && !(use.n <- all(is.n)))      use.n <- any(is.n <- knd == "n")
136        if(use.n && !(use.n <- all(is.n))) {
137          Tlst[is.n] <- lapply(Tlst[is.n], as, "lMatrix")          Tlst[is.n] <- lapply(Tlst[is.n], as, "lMatrix")
138            knd [is.n] <- "l"
139        }
140        use.l <- !use.n && all(knd == "l")
141      if(all(sym)) { ## result should be *symmetric*      if(all(sym)) { ## result should be *symmetric*
142          uplos <- sapply(Tlst, slot, "uplo") ## either "U" or "L"          uplos <- vapply(Tlst, slot, ".", "uplo") ## either "U" or "L"
143          tLU <- table(uplos)# of length 1 or 2 ..          tLU <- table(uplos)# of length 1 or 2 ..
144          if(length(tLU) == 1) { ## all "U" or all "L"          if(length(tLU) == 1) { ## all "U" or all "L"
145              useU <- uplos[1] == "U"              useU <- uplos[1] == "U"
# Line 121  Line 153 
153          if(use.n) { ## return nsparseMatrix :          if(use.n) { ## return nsparseMatrix :
154              r <- new("nsTMatrix")              r <- new("nsTMatrix")
155          } else {          } else {
156              r <- new("dsTMatrix")              r <- new(paste0(if(use.l) "l" else "d", "sTMatrix"))
157              r@x <- unlist(lapply(Tlst, slot, "x"))              r@x <- unlist(lapply(Tlst, slot, "x"))
158          }          }
159          r@uplo <- if(useU) "U" else "L"          r@uplo <- if(useU) "U" else "L"
160      }      }
161      else if(all(tri) && { ULs <- sapply(Tlst, slot, "uplo")##  "U" or "L"      else if(all(tri) && { ULs <- vapply(Tlst, slot, ".", "uplo")##  "U" or "L"
162                            all(ULs[1L] == ULs[-1L]) } ## all upper or all lower                            all(ULs[1L] == ULs[-1L]) } ## all upper or all lower
163         ){ ## *triangular* result         ){ ## *triangular* result
164    
165          if(use.n) { ## return nsparseMatrix :          if(use.n) { ## return nsparseMatrix :
166              r <- new("ntTMatrix")              r <- new("ntTMatrix")
167          } else {          } else {
168              r <- new("dtTMatrix")              r <- new(paste0(if(use.l) "l" else "d", "tTMatrix"))
169              r@x <- unlist(lapply(Tlst, slot, "x"))              r@x <- unlist(lapply(Tlst, slot, "x"))
170          }          }
171          r@uplo <- ULs[1L]          r@uplo <- ULs[1L]
# Line 144  Line 176 
176          if(use.n) { ## return nsparseMatrix :          if(use.n) { ## return nsparseMatrix :
177              r <- new("ngTMatrix")              r <- new("ngTMatrix")
178          } else {          } else {
179              r <- new("dgTMatrix")              r <- new(paste0(if(use.l) "l" else "d", "gTMatrix"))
180              r@x <- unlist(lapply(Tlst, slot, "x"))              r@x <- unlist(lapply(Tlst, slot, "x"))
181          }          }
182      }      }
# Line 170  Line 202 
202  .diag2tT <- function(from, uplo = "U", kind = .M.kind(from)) {  .diag2tT <- function(from, uplo = "U", kind = .M.kind(from)) {
203      ## to triangular Tsparse      ## to triangular Tsparse
204      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
205      new(paste(kind, "tTMatrix", sep=''),      new(paste0(kind, "tTMatrix"),
206          diag = from@diag, Dim = from@Dim, Dimnames = from@Dimnames,          diag = from@diag, Dim = from@Dim, Dimnames = from@Dimnames,
207          uplo = uplo,          uplo = uplo,
208          x = from@x, # <- ok for diag = "U" and "N" (!)          x = from@x, # <- ok for diag = "U" and "N" (!)
# Line 181  Line 213 
213      ## to symmetric Tsparse      ## to symmetric Tsparse
214      n <- from@Dim[1]      n <- from@Dim[1]
215      i <- seq_len(n) - 1L      i <- seq_len(n) - 1L
216      new(paste(kind, "sTMatrix", sep=''),      new(paste0(kind, "sTMatrix"),
217          Dim = from@Dim, Dimnames = from@Dimnames,          Dim = from@Dim, Dimnames = from@Dimnames,
218          i = i, j = i, uplo = uplo,          i = i, j = i, uplo = uplo,
219          x = if(from@diag == "N") from@x else ## "U"-diag          x = if(from@diag == "N") from@x else ## "U"-diag
# Line 190  Line 222 
222                         "l" =,                         "l" =,
223                         "n" = TRUE,                         "n" = TRUE,
224                         ## otherwise                         ## otherwise
225                         stop("'", kind,"' kind not yet implemented")), n))                         stop(gettextf("%s kind not yet implemented",
226                                         sQuote(kind)), domain=NA)),
227                    n))
228  }  }
229    
230  ## diagonal -> triangular,  upper / lower depending on "partner":  ## diagonal -> triangular,  upper / lower depending on "partner":
# Line 206  Line 240 
240          .diag2tT(d, uplo = if(extends(clx,"triangularMatrix")) x@uplo else "U", kind)          .diag2tT(d, uplo = if(extends(clx,"triangularMatrix")) x@uplo else "U", kind)
241  }  }
242    
243    ## FIXME: should not be needed {when ddi* is dsparse* etc}:
244    setMethod("is.finite", signature(x = "diagonalMatrix"),
245              function(x) is.finite(.diag2tT(x)))
246    setMethod("is.infinite", signature(x = "diagonalMatrix"),
247              function(x) is.infinite(.diag2tT(x)))
248    
249  ## In order to evade method dispatch ambiguity warnings,  ## In order to evade method dispatch ambiguity warnings,
250  ## 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
251  ## "hack"  instead of signature  x = "diagonalMatrix" :  ## "hack"  instead of signature  x = "diagonalMatrix" :
252  ##  ##
253  ## ddi*:  ## ddi*:
254  diag2tT <- function(from) .diag2tT(from, "U", "d")  di2tT <- function(from) .diag2tT(from, "U", "d")
255  setAs("ddiMatrix", "triangularMatrix", diag2tT)  setAs("ddiMatrix", "triangularMatrix", di2tT)
256  setAs("ddiMatrix", "sparseMatrix", diag2tT)  ##_no_longer_ setAs("ddiMatrix", "sparseMatrix", di2tT)
257  ## needed too (otherwise <dense> -> Tsparse is taken):  ## needed too (otherwise <dense> -> Tsparse is taken):
258  setAs("ddiMatrix", "TsparseMatrix", diag2tT)  setAs("ddiMatrix", "TsparseMatrix", di2tT)
259    setAs("ddiMatrix", "dsparseMatrix", di2tT)
260  setAs("ddiMatrix", "CsparseMatrix",  setAs("ddiMatrix", "CsparseMatrix",
261        function(from) as(.diag2tT(from, "U", "d"), "CsparseMatrix"))        function(from) as(.diag2tT(from, "U", "d"), "CsparseMatrix"))
262  setAs("ddiMatrix", "symmetricMatrix",  setAs("ddiMatrix", "symmetricMatrix",
263        function(from) .diag2sT(from, "U", "d"))        function(from) .diag2sT(from, "U", "d"))
264  ##  ##
265  ## ldi*:  ## ldi*:
266  diag2tT <- function(from) .diag2tT(from, "U", "l")  di2tT <- function(from) .diag2tT(from, "U", "l")
267  setAs("ldiMatrix", "triangularMatrix", diag2tT)  setAs("ldiMatrix", "triangularMatrix", di2tT)
268  setAs("ldiMatrix", "sparseMatrix", diag2tT)  ##_no_longer_ setAs("ldiMatrix", "sparseMatrix", di2tT)
269  ## needed too (otherwise <dense> -> Tsparse is taken):  ## needed too (otherwise <dense> -> Tsparse is taken):
270  setAs("ldiMatrix", "TsparseMatrix", diag2tT)  setAs("ldiMatrix", "TsparseMatrix", di2tT)
271    setAs("ldiMatrix", "lsparseMatrix", di2tT)
272  setAs("ldiMatrix", "CsparseMatrix",  setAs("ldiMatrix", "CsparseMatrix",
273        function(from) as(.diag2tT(from, "U", "l"), "CsparseMatrix"))        function(from) as(.diag2tT(from, "U", "l"), "CsparseMatrix"))
274  setAs("ldiMatrix", "symmetricMatrix",  setAs("ldiMatrix", "symmetricMatrix",
275        function(from) .diag2sT(from, "U", "l"))        function(from) .diag2sT(from, "U", "l"))
276    rm(di2tT)
277    
278  setAs("diagonalMatrix", "nMatrix",  setAs("diagonalMatrix", "nMatrix",
279        function(from) {        function(from) {
# Line 247  Line 288 
288  ## Cheap fast substitute for diag() which *does* preserve the mode of x :  ## Cheap fast substitute for diag() which *does* preserve the mode of x :
289  mkDiag <- function(x, n) {  mkDiag <- function(x, n) {
290      y <- matrix(as0(mod=mode(x)), n,n)      y <- matrix(as0(mod=mode(x)), n,n)
291      if (n > 0) y[1L + 0:(n - 1L) * (n + 1L)] <- x      if (n > 0) y[1L + 0:(n - 1L) * (n + 1)] <- x
292      y      y
293  }  }
294    
# Line 264  Line 305 
305                mod.x <- mode(x@x)                mod.x <- mode(x@x)
306                r <- vector(mod.x, length = n^2)                r <- vector(mod.x, length = n^2)
307                if(n)                if(n)
308                    r[1 + 0:(n - 1) * (n + 1)] <-                    r[1 + 0:(n - 1L) * (n + 1)] <-
309                        if(x@diag == "U") as1(mod=mod.x) else x@x                        if(x@diag == "U") as1(mod=mod.x) else x@x
310                r                r
311            })            })
# Line 272  Line 313 
313  setAs("diagonalMatrix", "generalMatrix", # prefer sparse:  setAs("diagonalMatrix", "generalMatrix", # prefer sparse:
314        function(from) as(as(from, "CsparseMatrix"), "generalMatrix"))        function(from) as(as(from, "CsparseMatrix"), "generalMatrix"))
315    
316    setAs("diagonalMatrix", "denseMatrix",
317          function(from) as(as(from, "CsparseMatrix"), "denseMatrix"))
318    
319    ..diag.x <- function(m)                   rep.int(as1(m@x), m@Dim[1])
320  .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
321    
322  .diag.2N <- function(m) {  .diag.2N <- function(m) {
# Line 292  Line 337 
337            d <- dim(from)            d <- dim(from)
338            if(d[1] != (n <- d[2])) stop("non-square matrix")            if(d[1] != (n <- d[2])) stop("non-square matrix")
339            if(any(from[row(from) != col(from)] != 0))            if(any(from[row(from) != col(from)] != 0))
340                stop("has non-zero off-diagonal entries")                stop("matrix with non-zero off-diagonals cannot be coerced to \"diagonalMatrix\"")
341            x <- diag(from)            x <- diag(from)
342            if(is.logical(x)) {            if(is.logical(x)) {
343                cl <- "ldiMatrix"                cl <- "ldiMatrix"
344                uni <- all(x)                uni <- allTrue(x) ## uni := {is it unit-diagonal ?}
345            } else {            } else {
346                cl <- "ddiMatrix"                cl <- "ddiMatrix"
347                uni <- all(x == 1)                uni <- allTrue(x == 1)
348                storage.mode(x) <- "double"                storage.mode(x) <- "double"
349            } ## TODO: complex            } ## TODO: complex
350            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 316  Line 361 
361            x <- diag(from)            x <- diag(from)
362            if(is.logical(x)) {            if(is.logical(x)) {
363                cl <- "ldiMatrix"                cl <- "ldiMatrix"
364                uni <- all(x)                uni <- allTrue(x)
365            } else {            } else {
366                cl <- "ddiMatrix"                cl <- "ddiMatrix"
367                uni <- all(x == 1)                uni <- allTrue(x == 1)
368                storage.mode(x) <- "double"                storage.mode(x) <- "double"
369            }            } ## TODO: complex
370            new(cl, Dim = c(n,n), diag = if(uni) "U" else "N",            new(cl, Dim = c(n,n), diag = if(uni) "U" else "N",
371                x = if(uni) x[FALSE] else x)                x = if(uni) x[FALSE] else x)
372        })        })
# Line 331  Line 376 
376            function(x = 1, nrow, ncol) .diag.x(x))            function(x = 1, nrow, ncol) .diag.x(x))
377    
378  subDiag <- function(x, i, j, ..., drop) {  subDiag <- function(x, i, j, ..., drop) {
379      x <- as(x, "sparseMatrix")      x <- as(x, "CsparseMatrix") ## << was "TsparseMatrix" (Csparse is faster now)
380      x <- if(missing(i))      x <- if(missing(i))
381          x[, j, drop=drop]          x[, j, drop=drop]
382      else if(missing(j))      else if(missing(j))
383          x[i, , drop=drop]          if(nargs() == 4) x[i, , drop=drop] else x[i, drop=drop]
384      else      else
385          x[i,j, drop=drop]          x[i,j, drop=drop]
386      if(isS4(x) && isDiagonal(x)) as(x, "diagonalMatrix") else x      if(isS4(x) && isDiagonal(x)) as(x, "diagonalMatrix") else x
# Line 345  Line 390 
390                           j = "index", drop = "logical"), subDiag)                           j = "index", drop = "logical"), subDiag)
391  setMethod("[", signature(x = "diagonalMatrix", i = "index",  setMethod("[", signature(x = "diagonalMatrix", i = "index",
392                          j = "missing", drop = "logical"),                          j = "missing", drop = "logical"),
393            function(x, i, j, ..., drop) subDiag(x, i=i, drop=drop))            function(x, i, j, ..., drop) {
394                  na <- nargs()
395                  Matrix.msg("diag[i,m,l] : nargs()=", na, .M.level = 2)
396                  if(na == 4)
397                       subDiag(x, i=i, , drop=drop)
398                  else subDiag(x, i=i,   drop=drop)
399              })
400  setMethod("[", signature(x = "diagonalMatrix", i = "missing",  setMethod("[", signature(x = "diagonalMatrix", i = "missing",
401                           j = "index", drop = "logical"),                           j = "index", drop = "logical"),
402            function(x, i, j, ..., drop) subDiag(x, j=j, drop=drop))            function(x, i, j, ..., drop) subDiag(x, j=j, drop=drop))
# Line 355  Line 406 
406  ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch  ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch
407  ## Only(?) current bug:  x[i] <- value  is wrong when  i is *vector*  ## Only(?) current bug:  x[i] <- value  is wrong when  i is *vector*
408  replDiag <- function(x, i, j, ..., value) {  replDiag <- function(x, i, j, ..., value) {
409      x <- as(x, "sparseMatrix")      x <- as(x, "CsparseMatrix")# was "Tsparse.." till 2012-07
410      if(missing(i))      if(missing(i))
411          x[, j] <- value          x[, j] <- value
412      else if(missing(j)) { ##  x[i , ] <- v  *OR*   x[i] <- v      else if(missing(j)) { ##  x[i , ] <- v  *OR*   x[i] <- v
# Line 365  Line 416 
416              x[i, ] <- value              x[i, ] <- value
417          else if(na == 3)          else if(na == 3)
418              x[i] <- value              x[i] <- value
419          else stop("Internal bug: nargs()=",na,"; please report")          else stop(gettextf("Internal bug: nargs()=%d; please report",
420                               na), domain=NA)
421      } else      } else
422          x[i,j] <- value          x[i,j] <- value
423      if(isDiagonal(x)) as(x, "diagonalMatrix") else x      if(isDiagonal(x)) as(x, "diagonalMatrix") else x
# Line 384  Line 436 
436                           replDiag(x, i=i, , value=value)                           replDiag(x, i=i, , value=value)
437                   })                   })
438    
439    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing",
440                                    j = "index", value = "replValue"),
441                     function(x,i,j, ..., value) replDiag(x, j=j, value=value))
442    
443    ## x[] <- value :
444    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing",
445                                    j = "missing", value = "ANY"),
446                     function(x,i,j, ..., value)
447                 {
448                  if(all0(value)) { # be faster
449                      r <- new(paste0(.M.kindC(getClassDef(class(x))),"tTMatrix"))# of all "0"
450                      r@Dim <- x@Dim
451                      r@Dimnames <- x@Dimnames
452                      r
453                  } else { ## typically non-sense: assigning to full sparseMatrix
454                      x[TRUE] <- value
455                      x
456                  }
457              })
458    
459    
460  setReplaceMethod("[", signature(x = "diagonalMatrix",  setReplaceMethod("[", signature(x = "diagonalMatrix",
461                                  i = "matrix", # 2-col.matrix                                  i = "matrix", # 2-col.matrix
462                                  j = "missing", value = "replValue"),                                  j = "missing", value = "replValue"),
# Line 395  Line 468 
468                                   if(any(value != one | is.na(value))) {                                   if(any(value != one | is.na(value))) {
469                                       x@diag <- "N"                                       x@diag <- "N"
470                                       x@x <- rep.int(one, x@Dim[1])                                       x@x <- rep.int(one, x@Dim[1])
471                                   }                                   } else return(x)
472                               }                               }
473                               x@x[ii] <- value                               x@x[ii] <- value
474                               x                               x
475                           } else { ## no longer diagonal, but remain sparse:                           } else { ## no longer diagonal, but remain sparse:
476                               x <- as(x, "sparseMatrix")                               x <- as(x, "TsparseMatrix")
477                               x[i] <- value                               x[i] <- value
478                               x                               x
479                           }                           }
# Line 412  Line 485 
485                       }                       }
486                   })                   })
487    
 setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing",  
                                 j = "index", value = "replValue"),  
                  function(x,i,j, ..., value) replDiag(x, j=j, value=value))  
488    
489    ## value = "sparseMatrix":
490  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",
491                                  value = "scarceMatrix"),                                  value = "sparseMatrix"),
492                   function (x, i, j, ..., value)                   function (x, i, j, ..., value)
493                   callGeneric(x=x, , j=j, value = as(value, "sparseVector")))                   callGeneric(x=x, , j=j, value = as(value, "sparseVector")))
494  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "missing",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "missing",
495                                  value = "scarceMatrix"),                                  value = "sparseMatrix"),
496                   function (x, i, j, ..., value)                   function (x, i, j, ..., value)
497                   callGeneric(x=x, i=i, , value = as(value, "sparseVector")))                   callGeneric(x=x, i=i, , value = as(value, "sparseVector")))
498  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "index",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "index",
499                                  value = "scarceMatrix"),                                  value = "sparseMatrix"),
500                   function (x, i, j, ..., value)                   function (x, i, j, ..., value)
501                   callGeneric(x=x, i=i, j=j, value = as(value, "sparseVector")))                   callGeneric(x=x, i=i, j=j, value = as(value, "sparseVector")))
502    
503    ## value = "sparseVector":
504  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",
505                                  value = "sparseVector"),                                  value = "sparseVector"),
506                   replDiag)                   replDiag)
# Line 443  Line 515 
515  setMethod("t", signature(x = "diagonalMatrix"),  setMethod("t", signature(x = "diagonalMatrix"),
516            function(x) { x@Dimnames <- x@Dimnames[2:1] ; x })            function(x) { x@Dimnames <- x@Dimnames[2:1] ; x })
517    
518  setMethod("isDiagonal", signature(object = "diagonalMatrix"),  setMethod("isDiagonal",   "diagonalMatrix", function(object) TRUE)
519            function(object) TRUE)  setMethod("isTriangular", "diagonalMatrix", function(object, ...) TRUE)
520  setMethod("isTriangular", signature(object = "diagonalMatrix"),  setMethod("isSymmetric",  "diagonalMatrix", function(object, ...) TRUE)
           function(object) TRUE)  
 setMethod("isSymmetric", signature(object = "diagonalMatrix"),  
           function(object, ...) TRUE)  
521    
522  setMethod("symmpart", signature(x = "diagonalMatrix"), function(x) x)  setMethod("symmpart", signature(x = "diagonalMatrix"), function(x) x)
523  setMethod("skewpart", signature(x = "diagonalMatrix"), setZero)  setMethod("skewpart", signature(x = "diagonalMatrix"), setZero)
# Line 485  Line 554 
554  ##       ---------------------  ##       ---------------------
555  ## Note that "ldi" logical are treated as numeric  ## Note that "ldi" logical are treated as numeric
556  diagdiagprod <- function(x, y) {  diagdiagprod <- function(x, y) {
557      n <- dimCheck(x,y)[1]      dimCheck(x,y)
558      if(x@diag != "U") {      if(x@diag != "U") {
559          if(y@diag != "U") {          if(y@diag != "U") {
560              nx <- x@x * y@x              nx <- x@x * y@x
# Line 517  Line 586 
586      dx <- dim(x)      dx <- dim(x)
587      dy <- dim(y)      dy <- dim(y)
588      if(dx[2] != dy[1]) stop("non-matching dimensions")      if(dx[2] != dy[1]) stop("non-matching dimensions")
     n <- dx[1]  
589      as(if(x@diag == "U") y else x@x * y, "Matrix")      as(if(x@diag == "U") y else x@x * y, "Matrix")
590  }  }
591  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),
# Line 586  Line 654 
654  ##        function(x, y = NULL) {  ##        function(x, y = NULL) {
655  ##           })  ##           })
656    
657  ## setMethod("tcrossprod", signature(x = "denseMatrix", y = "diagonalMatrix"),  ##' @param x CsparseMatrix
658  ##        function(x, y = NULL) {  ##' @param y diagonalMatrix
659  ##           })  ##' @return x %*% y
660    Cspdiagprod <- function(x, y) {
661        dx <- dim(x <- .Call(Csparse_diagU2N, x))
662        dy <- dim(y)
663        if(dx[2] != dy[1]) stop("non-matching dimensions")
664        if(y@diag == "N") { ## otherwise: y == Diagonal(n) : multiplication is identity
665            if(!all(y@x[1L] == y@x[-1L]) && is(x, "symmetricMatrix"))
666                x <- as(x, "generalMatrix")
667            ind <- rep.int(seq_len(dx[2]), x@p[-1] - x@p[-dx[2]-1L])
668            x@x <- x@x * y@x[ind]
669            if(.hasSlot(x, "factors") && length(x@factors)) {# drop cashed ones
670                ## instead of dropping all factors, be smart about some
671                ## TODO ......
672                x@factors <- list()
673            }
674        }
675        x
676    }
677    
678    ##' @param x diagonalMatrix
679    ##' @param y CsparseMatrix
680    ##' @return x %*% y
681    diagCspprod <- function(x, y) {
682        dx <- dim(x)
683        dy <- dim(y <- .Call(Csparse_diagU2N, y))
684        if(dx[2] != dy[1]) stop("non-matching dimensions")
685        if(x@diag == "N") {
686            if(!all(x@x[1L] == x@x[-1L]) && is(y, "symmetricMatrix"))
687                y <- as(y, "generalMatrix")
688            y@x <- y@x * x@x[y@i + 1L]
689            if(.hasSlot(y, "factors") && length(yf <- y@factors)) {
690                ## TODO
691                if(FALSE) { ## instead of dropping all factors, be smart about some
692                    keep <- character()
693                    if(any(iLU <- names(yf) == "LU")) {
694                        keep <- "LU"
695                    }
696                    y@factors <- yf[keep]
697                } else y@factors <- list() ## for now
698            }
699        }
700        y
701    }
702    
703    setMethod("crossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
704              function(x, y = NULL) diagCspprod(x, y))
705    
706  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
707            function(x, y = NULL) crossprod(as(x, "sparseMatrix"), y))            function(x, y = NULL) diagCspprod(x, as(y, "CsparseMatrix")))
708    
709    ## Prefer calling diagCspprod to Cspdiagprod if going to transpose anyway
710    ##  x'y == (y'x)'
711    setMethod("crossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
712              function(x, y = NULL) t(diagCspprod(y, x)))
713    
714  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
715            function(x, y = NULL) crossprod(x, as(y, "sparseMatrix")))            function(x, y = NULL) t(diagCspprod(y, as(x, "Csparsematrix"))))
716    
717    setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
718              function(x, y = NULL) diagCspprod(x, t(y)))
719    
720  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
721            function(x, y = NULL) tcrossprod(as(x, "sparseMatrix"), y))            function(x, y = NULL) diagCspprod(x, t(as(y, "CsparseMatrix"))))
722    
723    setMethod("tcrossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
724              function(x, y = NULL) Cspdiagprod(x, y))
725    
726  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
727            function(x, y = NULL) tcrossprod(x, as(y, "sparseMatrix")))            function(x, y = NULL) Cspdiagprod(as(x, "CsparseMatrix"), y))
728    
729    setMethod("%*%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
730              function(x, y) diagCspprod(x, y))
731    
732    ## instead of "sparseMatrix", use: [RT]sparse.. ("closer" in method dispatch)
733    for(cl in c("TsparseMatrix", "RsparseMatrix")) {
734    
 ## FIXME?: In theory, this can be done *FASTER*, in some cases, via tapply1()  
735  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
736            function(x, y) as(x, "sparseMatrix") %*% y)            function(x, y) diagCspprod(as(x, "CsparseMatrix"), y))
737    
738  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
739            function(x, y) x %*% as(y, "sparseMatrix"))            function(x, y) Cspdiagprod(as(x, "CsparseMatrix"), y))
740  ## NB: The previous is *not* triggering for  "ddi" o "dgC" (= distance 3)  }
741  ##     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)  
742  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
743            function(x, y) x %*% as(y, "CsparseMatrix"))            function(x, y) Cspdiagprod(x, y))
744  ## NB: this is *not* needed for Tsparse & Rsparse  
745  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*
746  ##       do indeed work by going through sparse (and *not* ddense)!  ##       do indeed work by going through sparse (and *not* ddense)!
747    
   
   
748  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),
749            function(a, b, ...) {            function(a, b, ...) {
750                a@x <- 1/ a@x                a@x <- 1/ a@x
# Line 629  Line 753 
753            })            })
754    
755  solveDiag <- function(a, b, ...) {  solveDiag <- function(a, b, ...) {
756      if((n <- a@Dim[1]) != nrow(b))      if(a@Dim[1] != nrow(b))
757          stop("incompatible matrix dimensions")          stop("incompatible matrix dimensions")
758      ## trivially invert a 'in place' and multiply:      ## trivially invert a 'in place' and multiply:
759      a@x <- 1/ a@x      a@x <- 1/ a@x
# Line 648  Line 772 
772  ###---------------- <Ops> (<Arith>, <Logic>, <Compare> ) ----------------------  ###---------------- <Ops> (<Arith>, <Logic>, <Compare> ) ----------------------
773    
774  ## 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.)  
775  diagOdiag <- function(e1,e2) {  diagOdiag <- function(e1,e2) {
776      ## result should also be diagonal _ if possible _      ## result should also be diagonal _ if possible _
777      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 656  Line 779 
779      r00 <- callGeneric(if(is.numeric(e1@x)) 0 else FALSE,      r00 <- callGeneric(if(is.numeric(e1@x)) 0 else FALSE,
780                         if(is.numeric(e2@x)) 0 else FALSE)                         if(is.numeric(e2@x)) 0 else FALSE)
781      if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* diagonal      if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* diagonal
782          if(is.numeric(r)) {          if(is.numeric(r)) { # "double" *or* "integer"
783              if(is.numeric(e2@x)) {              if(is.numeric(e2@x)) {
784                  e2@x <- r; return(.diag.2N(e2)) }                  e2@x <- r; return(.diag.2N(e2)) }
785              if(!is.numeric(e1@x))              if(!is.numeric(e1@x))
786                  ## e.g. e1, e2 are logical;                  ## e.g. e1, e2 are logical;
787                  e1 <- as(e1, "dMatrix")                  e1 <- as(e1, "dMatrix")
788                if(!is.double(r)) r <- as.double(r)
789          }          }
790          else if(is.logical(r))          else if(is.logical(r))
791              e1 <- as(e1, "lMatrix")              e1 <- as(e1, "lMatrix")
792          else stop("intermediate 'r' is of type", typeof(r))          else stop(gettextf("intermediate 'r' is of type %s",
793                               typeof(r)), domain=NA)
794          e1@x <- r          e1@x <- r
795          .diag.2N(e1)          .diag.2N(e1)
796      }      }
797      else { ## result not diagonal, but at least symmetric:      else { ## result not diagonal, but at least symmetric:
798            ## e.g., m == m
799          isNum <- (is.numeric(r) || is.numeric(r00))          isNum <- (is.numeric(r) || is.numeric(r00))
800          isLog <- (is.logical(r) || is.logical(r00))          isLog <- (is.logical(r) || is.logical(r00))
801            Matrix.msg("exploding <diag> o <diag> into dense matrix", .M.level = 2)
         if(getOption("verbose"))  
             message("exploding  <diag>  o  <diag>  into dense matrix")  
802          d <- e1@Dim          d <- e1@Dim
803          n <- d[1]          n <- d[1]
804          stopifnot(length(r) == n)          stopifnot(length(r) == n)
805            if(isNum && !is.double(r)) r <- as.double(r)
806          xx <- as.vector(matrix(rbind(r, matrix(r00,n,n)), n,n))          xx <- as.vector(matrix(rbind(r, matrix(r00,n,n)), n,n))
807          newcl <-          newcl <-
808              paste(if(isNum) "d" else if(isLog) {              paste0(if(isNum) "d" else if(isLog) {
809                  if(!any(is.na(r)) && !any(is.na(r00))) "n" else "l"                  if(!any(is.na(r)) && !any(is.na(r00))) "n" else "l"
810              } else stop("not yet implemented .. please report")              } else stop("not yet implemented .. please report"), "syMatrix")
                   ,  
                   "syMatrix", sep='')  
811    
812          new(newcl, Dim = e1@Dim, Dimnames = e1@Dimnames, x = xx)          new(newcl, Dim = e1@Dim, Dimnames = e1@Dimnames, x = xx)
813      }      }
# Line 702  Line 825 
825              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)
826  }  }
827    
828  ## FIXME:    diagonal  o  triangular  |-->  triangular  ## diagonal  o  triangular  |-->  triangular
829  ## -----     diagonal  o  symmetric   |-->  symmetric  ## diagonal  o  symmetric   |-->  symmetric
830  ##    {also when other is sparse: do these "here" --  ##    {also when other is sparse: do these "here" --
831  ##     before conversion to sparse, since that loses "diagonality"}  ##     before conversion to sparse, since that loses "diagonality"}
832    diagOtri <- function(e1,e2) {
833        ## result must be triangular
834        r <- callGeneric(d1 <- .diag.x(e1), diag(e2)) # error if not "compatible"
835        ## Check what happens with non-diagonals, i.e. (0 o 0), (FALSE o 0), ...:
836        e1.0 <- if(.n1 <- is.numeric(d1   )) 0 else FALSE
837        r00 <- callGeneric(e1.0, if(.n2 <- is.numeric(e2[0])) 0 else FALSE)
838        if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* triangular
839            diag(e2) <- r
840            ## check what happens "in the triangle"
841            e2.2 <- if(.n2) 2 else TRUE
842            if(!callGeneric(e1.0, e2.2) == e2.2) { # values "in triangle" can change:
843                n <- dim(e2)[1L]
844                it <- indTri(n, upper = (e2@uplo == "U"))
845                e2[it] <- callGeneric(e1.0, e2[it])
846            }
847            e2
848        }
849        else { ## result not triangular ---> general
850            rr <- as(e2, "generalMatrix")
851            diag(rr) <- r
852            rr
853        }
854    }
855    
856    
857    setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "triangularMatrix"),
858              diagOtri)
859    ## For the reverse,  Ops == "Arith" | "Compare" | "Logic"
860    ##   'Arith'  :=  '"+"', '"-"', '"*"', '"^"', '"%%"', '"%/%"', '"/"'
861    setMethod("Arith", signature(e1 = "triangularMatrix", e2 = "diagonalMatrix"),
862              function(e1,e2)
863          { ## this must only trigger for *dense* e1
864              switch(.Generic,
865                     "+" = .Call(dtrMatrix_addDiag, as(e1,"dtrMatrix"),   .diag.x(e2)),
866                     "-" = .Call(dtrMatrix_addDiag, as(e1,"dtrMatrix"), - .diag.x(e2)),
867                     "*" = {
868                         n <- e2@Dim[1L]
869                         d2 <- if(e2@diag == "U") { # unit-diagonal
870                             d <- rep.int(as1(e2@x), n)
871                             e2@x <- d
872                             e2@diag <- "N"
873                             d
874                         } else e2@x
875                         e2@x <- diag(e1) * d2
876                         e2
877                     },
878                     "^" = { ## will be dense ( as  <ANY> ^ 0 == 1 ):
879                         e1 ^ as(e2, "denseMatrix")
880                     },
881                     ## otherwise:
882                     callGeneric(e1, diag2Tsmart(e2,e1)))
883    })
884    
885    ## Compare --> 'swap' (e.g.   e1 < e2   <==>  e2 > e1 ):
886    setMethod("Compare", signature(e1 = "triangularMatrix", e2 = "diagonalMatrix"),
887              .Cmp.swap)
888    ## '&' and "|'  are commutative:
889    setMethod("Logic", signature(e1 = "triangularMatrix", e2 = "diagonalMatrix"),
890              function(e1,e2) callGeneric(e2, e1))
891    
892  ## For almost everything else, diag* shall be treated "as sparse" :  ## For almost everything else, diag* shall be treated "as sparse" :
893  ## These are cheap implementations via coercion  ## These are cheap implementations via coercion
# Line 729  Line 911 
911  ##       Compare --> logical  ##       Compare --> logical
912  ##       Logic   --> logical: "lMatrix"  ##       Logic   --> logical: "lMatrix"
913    
914  ##  other = "numeric" : stay diagonal if possible  ## Other = "numeric" : stay diagonal if possible
915  ## ddi*: Arith: result numeric, potentially ddiMatrix  ## ddi*: Arith: result numeric, potentially ddiMatrix
916  setMethod("Arith", signature(e1 = "ddiMatrix", e2 = "numeric"),  for(arg2 in c("numeric","logical"))
917    setMethod("Arith", signature(e1 = "ddiMatrix", e2 = arg2),
918            function(e1,e2) {            function(e1,e2) {
919                n <- e1@Dim[1]; nsq <- n^2                n <- e1@Dim[1]
920                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
921                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
922                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
923                    if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)                    if(e1@diag == "U") {
924                    if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {                        if(any((r <- callGeneric(1, e2)) != 1)) {
925                        e1@diag <- "N"                        e1@diag <- "N"
926                        if(L1) r <- rep.int(r, n)                            e1@x[seq_len(n)] <- r # possibly recycling r
927                    } else                        } ## else: result = e1  (is "U" diag)
928                      } else {
929                        r <- callGeneric(e1@x, e2)                        r <- callGeneric(e1@x, e2)
930                    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
931                    return(e1)                        e1@x[] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
932                }                }
933                      e1
934                  } else
935                callGeneric(diag2tT.u(e1,e2, "d"), e2)                callGeneric(diag2tT.u(e1,e2, "d"), e2)
936            })            })
937    
938  setMethod("Arith", signature(e1 = "numeric", e2 = "ddiMatrix"),  for(arg1 in c("numeric","logical"))
939    setMethod("Arith", signature(e1 = arg1, e2 = "ddiMatrix"),
940            function(e1,e2) {            function(e1,e2) {
941                n <- e2@Dim[1]; nsq <- n^2                n <- e2@Dim[1]
942                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
943                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
944                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
945                    if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)                    if(e2@diag == "U") {
946                    if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {                        if(any((r <- callGeneric(e1, 1)) != 1)) {
947                        e2@diag <- "N"                        e2@diag <- "N"
948                        if(L1) r <- rep.int(r, n)                            e2@x[seq_len(n)] <- r # possibly recycling r
949                    } else                        } ## else: result = e2  (is "U" diag)
950                      } else {
951                        r <- callGeneric(e1, e2@x)                        r <- callGeneric(e1, e2@x)
952                    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
953                    return(e2)                        e2@x[] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
954                }                }
955                      e2
956                  } else
957                callGeneric(e1, diag2tT.u(e2,e1, "d"))                callGeneric(e1, diag2tT.u(e2,e1, "d"))
958            })            })
959    
960  ## ldi* Arith --> result numeric, potentially ddiMatrix  ## ldi* Arith --> result numeric, potentially ddiMatrix
961  setMethod("Arith", signature(e1 = "ldiMatrix", e2 = "numeric"),  for(arg2 in c("numeric","logical"))
962    setMethod("Arith", signature(e1 = "ldiMatrix", e2 = arg2),
963            function(e1,e2) {            function(e1,e2) {
964                n <- e1@Dim[1]; nsq <- n^2                n <- e1@Dim[1]
965                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
966                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
967                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
968                    if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)                    E <- copyClass(e1, "ddiMatrix", c("diag", "Dim", "Dimnames"), check=FALSE)
969                    if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {                    ## storage.mode(E@x) <- "double"
970                        e1@diag <- "N"                    if(e1@diag == "U") {
971                        if(L1) r <- rep.int(r, n)                        if(any((r <- callGeneric(1, e2)) != 1)) {
972                    } else                            E@diag <- "N"
973                              E@x[seq_len(n)] <- r # possibly recycling r
974                          } ## else: result = E  (is "U" diag)
975                      } else {
976                        r <- callGeneric(e1@x, e2)                        r <- callGeneric(e1@x, e2)
977                    e1 <- copyClass(e1, "ddiMatrix", c("diag", "Dim", "Dimnames"))                        ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
978                    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)  
979                }                }
980                callGeneric(diag2tT.u(e1,e2, "d"), e2)                    E
981                  } else
982                      callGeneric(diag2tT.u(e1,e2, "l"), e2)
983            })            })
984    
985  setMethod("Arith", signature(e1 = "numeric", e2 = "ldiMatrix"),  for(arg1 in c("numeric","logical"))
986    setMethod("Arith", signature(e1 = arg1, e2 = "ldiMatrix"),
987            function(e1,e2) {            function(e1,e2) {
988                n <- e2@Dim[1]; nsq <- n^2                n <- e2@Dim[1]
989                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
990                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
991                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
992                    if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)                    E <- copyClass(e2, "ddiMatrix", c("diag", "Dim", "Dimnames"), check=FALSE)
993                    if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {                    ## storage.mode(E@x) <- "double"
994                        e2@diag <- "N"                    if(e2@diag == "U") {
995                        if(L1) r <- rep.int(r, n)                        if(any((r <- callGeneric(e1, 1)) != 1)) {
996                    } else                            E@diag <- "N"
997                              E@x[seq_len(n)] <- r # possibly recycling r
998                          } ## else: result = E  (is "U" diag)
999                      } else {
1000                        r <- callGeneric(e1, e2@x)                        r <- callGeneric(e1, e2@x)
1001                    e2 <- copyClass(e2, "ddiMatrix", c("diag", "Dim", "Dimnames"))                        ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
1002                    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)  
1003                }                }
1004                callGeneric(e1, diag2tT.u(e2,e1, "d"))                    E
1005                  } else
1006                      callGeneric(e1, diag2tT.u(e2,e1, "l"))
1007            })            })
1008    
1009  ## ddi*: for "Ops" without Arith --> result logical, potentially ldi  ## ddi*: for "Ops" without "Arith": <Compare> or <Logic> --> result logical, potentially ldi
1010  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "numeric"),  ##
1011    ## Note that  ("numeric", "ddiMatrix")  is simply swapped, e.g.,
1012    if(FALSE) {
1013        selectMethod("<", c("numeric","lMatrix"))# Compare
1014        selectMethod("&", c("numeric","lMatrix"))# Logic
1015    } ## so we don't need to define a method here :
1016    
1017    for(arg2 in c("numeric","logical"))
1018    setMethod("Ops", signature(e1 = "ddiMatrix", e2 = arg2),
1019            function(e1,e2) {            function(e1,e2) {
1020                n <- e1@Dim[1]; nsq <- n^2                n <- e1@Dim[1]
1021                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
1022                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
1023                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
1024                    if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)                    E <- copyClass(e1, "ldiMatrix", c("diag", "Dim", "Dimnames"), check=FALSE)
1025                    if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {                    ## storage.mode(E@x) <- "logical"
1026                        e1@diag <- "N"                    if(e1@diag == "U") {
1027                        if(L1) r <- rep.int(r, n)                        if(any((r <- callGeneric(1, e2)) != 1)) {
1028                    } else                            E@diag <- "N"
1029                              E@x[seq_len(n)] <- r # possibly recycling r
1030                          } ## else: result = E  (is "U" diag)
1031                      } else {
1032                        r <- callGeneric(e1@x, e2)                        r <- callGeneric(e1@x, e2)
1033                    e1 <- copyClass(e1, "ldiMatrix", c("diag", "Dim", "Dimnames"))                        ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
1034                    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)  
1035                }                }
1036                callGeneric(diag2tT.u(e1,e2, "l"), e2)                    E
           })  
   
 setMethod("Ops", signature(e1 = "numeric", e2 = "ddiMatrix"),  
           function(e1,e2) {  
               n <- e2@Dim[1]; nsq <- n^2  
               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)  
1037                    } else                    } else
1038                        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"))  
1039            })            })
1040    
1041  ## ldi*: for "Ops" without Arith --> result logical, potentially ldi  ## ldi*: for "Ops" without "Arith": <Compare> or <Logic> --> result logical, potentially ldi
1042  setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "numeric"),  for(arg2 in c("numeric","logical"))
1043    setMethod("Ops", signature(e1 = "ldiMatrix", e2 = arg2),
1044            function(e1,e2) {            function(e1,e2) {
1045                n <- e1@Dim[1]; nsq <- n^2                n <- e1@Dim[1]
1046                f0 <- callGeneric(FALSE, e2)                f0 <- callGeneric(FALSE, e2)
1047                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
1048                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
1049                    if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)  
1050                    if(e1@diag == "U" && any((r <- callGeneric(TRUE, e2)) != 1)) {                    if(e1@diag == "U") {
1051                          if(any((r <- callGeneric(TRUE, e2)) != 1)) {
1052                        e1@diag <- "N"                        e1@diag <- "N"
1053                        if(L1) r <- rep.int(r, n)                            e1@x[seq_len(n)] <- r # possibly recycling r
1054                    } else                        } ## else: result = e1  (is "U" diag)
1055                      } else {
1056                        r <- callGeneric(e1@x, e2)                        r <- callGeneric(e1@x, e2)
1057                    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
1058                    return(e1)                        e1@x[] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
1059                }                }
1060                callGeneric(diag2tT.u(e1,e2, "l"), e2)                    e1
           })  
   
 setMethod("Ops", signature(e1 = "numeric", e2 = "ldiMatrix"),  
           function(e1,e2) {  
               n <- e2@Dim[1]; nsq <- n^2  
               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)  
1061                    } else                    } else
1062                        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"))  
1063            })            })
1064    
1065    
   
1066  ## Not {"sparseMatrix", "numeric} :  {"denseMatrix", "matrix", ... }  ## Not {"sparseMatrix", "numeric} :  {"denseMatrix", "matrix", ... }
1067  for(other in c("ANY", "Matrix", "dMatrix")) {  for(other in c("ANY", "Matrix", "dMatrix")) {
1068      ## ddi*:      ## ddi*:
1069      setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other),      setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other),
1070                function(e1,e2) callGeneric(diag2tT.u(e1,e2, "d"), e2))                function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2))
1071      setMethod("Ops", signature(e1 = other, e2 = "ddiMatrix"),      setMethod("Ops", signature(e1 = other, e2 = "ddiMatrix"),
1072                function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "d")))                function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "d")))
1073      ## ldi*:      ## ldi*:
1074      setMethod("Ops", signature(e1 = "ldiMatrix", e2 = other),      setMethod("Ops", signature(e1 = "ldiMatrix", e2 = other),
1075                function(e1,e2) callGeneric(diag2tT.u(e1,e2, "l"), e2))                function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2))
1076      setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),      setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),
1077                function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "l")))                function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "l")))
1078  }  }
1079  ## This should *not* dispatch to <dense> methods (in ./Ops.R ), as  
1080  ##  FALSE & <anything> |-> FALSE : hence result should be diagonal:  ## Direct subclasses of "denseMatrix": currently ddenseMatrix, ldense... :
1081  for(cl in diCls) {  dense.subCl <- local({ dM.scl <- getClass("denseMatrix")@subclasses
1082      setMethod("&", signature(e1 = cl, e2 = "ANY"),                         names(dM.scl)[vapply(dM.scl, slot, 0, "distance") == 1] })
1083                function(e1,e2) e1 & as(e2,"Matrix"))  for(DI in diCls) {
1084      setMethod("&", signature(e1 = "ANY", e2 = cl),      dMeth <- if(extends(DI, "dMatrix"))
1085                function(e1,e2) as(e1,"Matrix") & e2)          function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2)
1086      for(c2 in c("denseMatrix", "Matrix")) {      else # "lMatrix", the only other kind for now
1087          setMethod("&", signature(e1 = cl, e2 = c2),          function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2)
1088                    function(e1,e2) e1 & Diagonal(x = diag(e2)))      for(c2 in c(dense.subCl, "Matrix")) {
1089          setMethod("&", signature(e1 = c2, e2 = cl),          for(Fun in c("*", "&")) {
1090                    function(e1,e2) Diagonal(x = diag(e1)) & e2)              setMethod(Fun, signature(e1 = DI, e2 = c2),
1091                          function(e1,e2) callGeneric(e1, Diagonal(x = diag(e2))))
1092                setMethod(Fun, signature(e1 = c2, e2 = DI),
1093                          function(e1,e2) callGeneric(Diagonal(x = diag(e1)), e2))
1094            }
1095            setMethod("^", signature(e1 = c2, e2 = DI),
1096                      function(e1,e2) callGeneric(Diagonal(x = diag(e1)), e2))
1097            for(Fun in c("%%", "%/%", "/")) ## 0 <op> 0 |--> NaN  for these.
1098                setMethod(Fun, signature(e1 = DI, e2 = c2), dMeth)
1099      }      }
1100  }  }
1101    
1102    
1103  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"
1104  ### ----------  any, all: separately here  ### ----------   the last 4: separately here
1105  for(cl in diCls) {  for(cl in diCls) {
1106  setMethod("any", cl,  setMethod("any", cl,
1107            function (x, ..., na.rm) {            function (x, ..., na.rm) {
1108                if(any(x@Dim == 0)) FALSE                if(any(x@Dim == 0)) FALSE
1109                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)
1110            })            })
1111  setMethod("all",  cl, function (x, ..., na.rm) any(x@Dim == 0))  setMethod("all",  cl, function (x, ..., na.rm) {
1112  setMethod("prod", cl, function (x, ..., na.rm) as.numeric(any(x@Dim == 0)))      n <- x@Dim[1]
1113        if(n >= 2) FALSE
1114        else if(n == 0 || x@diag == "U") TRUE
1115        else all(x@x, ..., na.rm = na.rm)
1116    })
1117    setMethod("prod", cl, function (x, ..., na.rm) {
1118        n <- x@Dim[1]
1119        if(n >= 2) 0
1120        else if(n == 0 || x@diag == "U") 1
1121        else ## n == 1, diag = "N" :
1122            prod(x@x, ..., na.rm = na.rm)
1123    })
1124    
1125  setMethod("sum", cl,  setMethod("sum", cl,
1126            function(x, ..., na.rm) {            function(x, ..., na.rm) {
# Line 979  Line 1177 
1177                    invisible(object)                    invisible(object)
1178                }                }
1179            })            })
1180    
1181    rm(arg1, arg2, other, DI, cl, c1, c2,
1182       dense.subCl, diCls)# not used elsewhere
1183    
1184    setMethod("summary", signature(object = "diagonalMatrix"),
1185              function(object, ...) {
1186                  d <- dim(object)
1187                  r <- summary(object@x, ...)
1188                  attr(r, "header") <-
1189                      sprintf('%d x %d diagonal Matrix of class "%s"',
1190                              d[1], d[2], class(object))
1191                  ## use ole' S3 technology for such a simple case
1192                  class(r) <- c("diagSummary", class(r))
1193                  r
1194              })
1195    
1196    print.diagSummary <- function (x, ...) {
1197        cat(attr(x, "header"),"\n")
1198        class(x) <- class(x)[-1]
1199        print(x, ...)
1200        invisible(x)
1201    }

Legend:
Removed from v.2237  
changed lines
  Added in v.2984

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