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 2891, Fri Aug 16 14:09:01 2013 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
# Line 214  Line 253 
253  ## ddi*:  ## ddi*:
254  diag2tT <- function(from) .diag2tT(from, "U", "d")  diag2tT <- function(from) .diag2tT(from, "U", "d")
255  setAs("ddiMatrix", "triangularMatrix", diag2tT)  setAs("ddiMatrix", "triangularMatrix", diag2tT)
256  setAs("ddiMatrix", "sparseMatrix", diag2tT)  ##_no_longer_ setAs("ddiMatrix", "sparseMatrix", diag2tT)
257  ## needed too (otherwise <dense> -> Tsparse is taken):  ## needed too (otherwise <dense> -> Tsparse is taken):
258  setAs("ddiMatrix", "TsparseMatrix", diag2tT)  setAs("ddiMatrix", "TsparseMatrix", diag2tT)
259    setAs("ddiMatrix", "dsparseMatrix", diag2tT)
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",
# Line 225  Line 265 
265  ## ldi*:  ## ldi*:
266  diag2tT <- function(from) .diag2tT(from, "U", "l")  diag2tT <- function(from) .diag2tT(from, "U", "l")
267  setAs("ldiMatrix", "triangularMatrix", diag2tT)  setAs("ldiMatrix", "triangularMatrix", diag2tT)
268  setAs("ldiMatrix", "sparseMatrix", diag2tT)  ##_no_longer_ setAs("ldiMatrix", "sparseMatrix", diag2tT)
269  ## needed too (otherwise <dense> -> Tsparse is taken):  ## needed too (otherwise <dense> -> Tsparse is taken):
270  setAs("ldiMatrix", "TsparseMatrix", diag2tT)  setAs("ldiMatrix", "TsparseMatrix", diag2tT)
271    setAs("ldiMatrix", "lsparseMatrix", diag2tT)
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",
# 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 485  Line 557 
557  ##       ---------------------  ##       ---------------------
558  ## Note that "ldi" logical are treated as numeric  ## Note that "ldi" logical are treated as numeric
559  diagdiagprod <- function(x, y) {  diagdiagprod <- function(x, y) {
560      n <- dimCheck(x,y)[1]      dimCheck(x,y)
561      if(x@diag != "U") {      if(x@diag != "U") {
562          if(y@diag != "U") {          if(y@diag != "U") {
563              nx <- x@x * y@x              nx <- x@x * y@x
# Line 517  Line 589 
589      dx <- dim(x)      dx <- dim(x)
590      dy <- dim(y)      dy <- dim(y)
591      if(dx[2] != dy[1]) stop("non-matching dimensions")      if(dx[2] != dy[1]) stop("non-matching dimensions")
     n <- dx[1]  
592      as(if(x@diag == "U") y else x@x * y, "Matrix")      as(if(x@diag == "U") y else x@x * y, "Matrix")
593  }  }
594  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),
# Line 586  Line 657 
657  ##        function(x, y = NULL) {  ##        function(x, y = NULL) {
658  ##           })  ##           })
659    
660  ## setMethod("tcrossprod", signature(x = "denseMatrix", y = "diagonalMatrix"),  Cspdiagprod <- function(x, y) {
661  ##        function(x, y = NULL) {      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") {
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        }
670        if(is(x, "compMatrix") && length(xf <- x@factors)) {
671            ## instead of dropping all factors, be smart about some
672            ## TODO ......
673            x@factors <- list()
674        }
675        x
676    }
677    
678    diagCspprod <- function(x, y) {
679        dx <- dim(x)
680        dy <- dim(y <- .Call(Csparse_diagU2N, y))
681        if(dx[2] != dy[1]) stop("non-matching dimensions")
682        if(x@diag == "N") {
683            if(!all(x@x[1L] == x@x[-1L]) && is(y, "symmetricMatrix"))
684                y <- as(y, "generalMatrix")
685            y@x <- y@x * x@x[y@i + 1L]
686        }
687        if(is(y, "compMatrix") && length(yf <- y@factors)) {
688            ## instead of dropping all factors, be smart about some
689            ## TODO
690            keep <- character()
691            if(iLU <- names(yf) == "LU") {
692                ## TODO keep <- "LU"
693            }
694            y@factors <- yf[keep]
695        }
696        y
697    }
698    
699    setMethod("crossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
700              function(x, y = NULL) diagCspprod(x, y))
701    
702  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
703            function(x, y = NULL) crossprod(as(x, "sparseMatrix"), y))            function(x, y = NULL) diagCspprod(x, as(y, "CsparseMatrix")))
704    
705    ## Prefer calling diagCspprod to Cspdiagprod if going to transpose anyway
706    ##  x'y == (y'x)'
707    setMethod("crossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
708              function(x, y = NULL) t(diagCspprod(y, x)))
709    
710  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
711            function(x, y = NULL) crossprod(x, as(y, "sparseMatrix")))            function(x, y = NULL) t(diagCspprod(y, as(x, "Csparsematrix"))))
712    
713    setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
714              function(x, y = NULL) diagCspprod(x, t(y)))
715    
716  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
717            function(x, y = NULL) tcrossprod(as(x, "sparseMatrix"), y))            function(x, y = NULL) diagCspprod(x, t(as(y, "CsparseMatrix"))))
718    
719    setMethod("tcrossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
720              function(x, y = NULL) Cspdiagprod(x, y))
721    
722  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
723            function(x, y = NULL) tcrossprod(x, as(y, "sparseMatrix")))            function(x, y = NULL) Cspdiagprod(as(x, "CsparseMatrix"), y))
724    
725    setMethod("%*%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
726              function(x, y) diagCspprod(x, y))
727    
 ## FIXME?: In theory, this can be done *FASTER*, in some cases, via tapply1()  
728  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
729            function(x, y) as(x, "sparseMatrix") %*% y)            function(x, y) diagCspprod(as(x, "CsparseMatrix"), y))
730    
731  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
732            function(x, y) x %*% as(y, "sparseMatrix"))            function(x, y) Cspdiagprod(as(x, "CsparseMatrix"), y))
733  ## 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)  
734  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
735            function(x, y) x %*% as(y, "CsparseMatrix"))            function(x, y) Cspdiagprod(x, y))
736  ## NB: this is *not* needed for Tsparse & Rsparse  
737  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*
738  ##       do indeed work by going through sparse (and *not* ddense)!  ##       do indeed work by going through sparse (and *not* ddense)!
739    
   
   
740  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),
741            function(a, b, ...) {            function(a, b, ...) {
742                a@x <- 1/ a@x                a@x <- 1/ a@x
# Line 629  Line 745 
745            })            })
746    
747  solveDiag <- function(a, b, ...) {  solveDiag <- function(a, b, ...) {
748      if((n <- a@Dim[1]) != nrow(b))      if(a@Dim[1] != nrow(b))
749          stop("incompatible matrix dimensions")          stop("incompatible matrix dimensions")
750      ## trivially invert a 'in place' and multiply:      ## trivially invert a 'in place' and multiply:
751      a@x <- 1/ a@x      a@x <- 1/ a@x
# Line 648  Line 764 
764  ###---------------- <Ops> (<Arith>, <Logic>, <Compare> ) ----------------------  ###---------------- <Ops> (<Arith>, <Logic>, <Compare> ) ----------------------
765    
766  ## 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.)  
767  diagOdiag <- function(e1,e2) {  diagOdiag <- function(e1,e2) {
768      ## result should also be diagonal _ if possible _      ## result should also be diagonal _ if possible _
769      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 771 
771      r00 <- callGeneric(if(is.numeric(e1@x)) 0 else FALSE,      r00 <- callGeneric(if(is.numeric(e1@x)) 0 else FALSE,
772                         if(is.numeric(e2@x)) 0 else FALSE)                         if(is.numeric(e2@x)) 0 else FALSE)
773      if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* diagonal      if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* diagonal
774          if(is.numeric(r)) {          if(is.numeric(r)) { # "double" *or* "integer"
775              if(is.numeric(e2@x)) {              if(is.numeric(e2@x)) {
776                  e2@x <- r; return(.diag.2N(e2)) }                  e2@x <- r; return(.diag.2N(e2)) }
777              if(!is.numeric(e1@x))              if(!is.numeric(e1@x))
778                  ## e.g. e1, e2 are logical;                  ## e.g. e1, e2 are logical;
779                  e1 <- as(e1, "dMatrix")                  e1 <- as(e1, "dMatrix")
780                if(!is.double(r)) r <- as.double(r)
781          }          }
782          else if(is.logical(r))          else if(is.logical(r))
783              e1 <- as(e1, "lMatrix")              e1 <- as(e1, "lMatrix")
784          else stop("intermediate 'r' is of type", typeof(r))          else stop(gettextf("intermediate 'r' is of type %s",
785                               typeof(r)), domain=NA)
786          e1@x <- r          e1@x <- r
787          .diag.2N(e1)          .diag.2N(e1)
788      }      }
789      else { ## result not diagonal, but at least symmetric:      else { ## result not diagonal, but at least symmetric:
790            ## e.g., m == m
791          isNum <- (is.numeric(r) || is.numeric(r00))          isNum <- (is.numeric(r) || is.numeric(r00))
792          isLog <- (is.logical(r) || is.logical(r00))          isLog <- (is.logical(r) || is.logical(r00))
793            Matrix.msg("exploding <diag> o <diag> into dense matrix", .M.level = 2)
         if(getOption("verbose"))  
             message("exploding  <diag>  o  <diag>  into dense matrix")  
794          d <- e1@Dim          d <- e1@Dim
795          n <- d[1]          n <- d[1]
796          stopifnot(length(r) == n)          stopifnot(length(r) == n)
797            if(isNum && !is.double(r)) r <- as.double(r)
798          xx <- as.vector(matrix(rbind(r, matrix(r00,n,n)), n,n))          xx <- as.vector(matrix(rbind(r, matrix(r00,n,n)), n,n))
799          newcl <-          newcl <-
800              paste(if(isNum) "d" else if(isLog) {              paste0(if(isNum) "d" else if(isLog) {
801                  if(!any(is.na(r)) && !any(is.na(r00))) "n" else "l"                  if(!any(is.na(r)) && !any(is.na(r00))) "n" else "l"
802              } else stop("not yet implemented .. please report")              } else stop("not yet implemented .. please report"), "syMatrix")
                   ,  
                   "syMatrix", sep='')  
803    
804          new(newcl, Dim = e1@Dim, Dimnames = e1@Dimnames, x = xx)          new(newcl, Dim = e1@Dim, Dimnames = e1@Dimnames, x = xx)
805      }      }
# Line 702  Line 817 
817              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)
818  }  }
819    
820  ## FIXME:    diagonal  o  triangular  |-->  triangular  ## diagonal  o  triangular  |-->  triangular
821  ## -----     diagonal  o  symmetric   |-->  symmetric  ## diagonal  o  symmetric   |-->  symmetric
822  ##    {also when other is sparse: do these "here" --  ##    {also when other is sparse: do these "here" --
823  ##     before conversion to sparse, since that loses "diagonality"}  ##     before conversion to sparse, since that loses "diagonality"}
824    diagOtri <- function(e1,e2) {
825        ## result must be triangular
826        r <- callGeneric(d1 <- .diag.x(e1), diag(e2)) # error if not "compatible"
827        ## Check what happens with non-diagonals, i.e. (0 o 0), (FALSE o 0), ...:
828        e1.0 <- if(.n1 <- is.numeric(d1   )) 0 else FALSE
829        r00 <- callGeneric(e1.0, if(.n2 <- is.numeric(e2[0])) 0 else FALSE)
830        if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* triangular
831            diag(e2) <- r
832            ## check what happens "in the triangle"
833            e2.2 <- if(.n2) 2 else TRUE
834            if(!callGeneric(e1.0, e2.2) == e2.2) { # values "in triangle" can change:
835                n <- dim(e2)[1L]
836                it <- indTri(n, upper = (e2@uplo == "U"))
837                e2[it] <- callGeneric(e1.0, e2[it])
838            }
839            e2
840        }
841        else { ## result not triangular ---> general
842            rr <- as(e2, "generalMatrix")
843            diag(rr) <- r
844            rr
845        }
846    }
847    
848    
849    setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "triangularMatrix"),
850              diagOtri)
851    ## For the reverse,  Ops == "Arith" | "Compare" | "Logic"
852    ##   'Arith'  :=  '"+"', '"-"', '"*"', '"^"', '"%%"', '"%/%"', '"/"'
853    setMethod("Arith", signature(e1 = "triangularMatrix", e2 = "diagonalMatrix"),
854              function(e1,e2)
855          { ## this must only trigger for *dense* e1
856              switch(.Generic,
857                     "+" = .Call(dtrMatrix_addDiag, as(e1,"dtrMatrix"),   .diag.x(e2)),
858                     "-" = .Call(dtrMatrix_addDiag, as(e1,"dtrMatrix"), - .diag.x(e2)),
859                     "*" = {
860                         n <- e2@Dim[1L]
861                         d2 <- if(e2@diag == "U") { # unit-diagonal
862                             d <- rep.int(as1(e2@x), n)
863                             e2@x <- d
864                             e2@diag <- "N"
865                             d
866                         } else e2@x
867                         e2@x <- diag(e1) * d2
868                         e2
869                     },
870                     "^" = { ## will be dense ( as  <ANY> ^ 0 == 1 ):
871                         e1 ^ as(e2, "denseMatrix")
872                     },
873                     ## otherwise:
874                     callGeneric(e1, diag2Tsmart(e2,e1)))
875    })
876    
877    ## Compare --> 'swap' (e.g.   e1 < e2   <==>  e2 > e1 ):
878    setMethod("Compare", signature(e1 = "triangularMatrix", e2 = "diagonalMatrix"),
879              .Cmp.swap)
880    ## '&' and "|'  are commutative:
881    setMethod("Logic", signature(e1 = "triangularMatrix", e2 = "diagonalMatrix"),
882              function(e1,e2) callGeneric(e2, e1))
883    
884  ## For almost everything else, diag* shall be treated "as sparse" :  ## For almost everything else, diag* shall be treated "as sparse" :
885  ## These are cheap implementations via coercion  ## These are cheap implementations via coercion
# Line 729  Line 903 
903  ##       Compare --> logical  ##       Compare --> logical
904  ##       Logic   --> logical: "lMatrix"  ##       Logic   --> logical: "lMatrix"
905    
906  ##  other = "numeric" : stay diagonal if possible  ## Other = "numeric" : stay diagonal if possible
907  ## ddi*: Arith: result numeric, potentially ddiMatrix  ## ddi*: Arith: result numeric, potentially ddiMatrix
908  setMethod("Arith", signature(e1 = "ddiMatrix", e2 = "numeric"),  for(arg2 in c("numeric","logical"))
909    setMethod("Arith", signature(e1 = "ddiMatrix", e2 = arg2),
910            function(e1,e2) {            function(e1,e2) {
911                n <- e1@Dim[1]; nsq <- n^2                n <- e1@Dim[1]
912                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
913                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
914                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
915                    if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)                    if(e1@diag == "U") {
916                    if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {                        if(any((r <- callGeneric(1, e2)) != 1)) {
917                        e1@diag <- "N"                        e1@diag <- "N"
918                        if(L1) r <- rep.int(r, n)                            e1@x[seq_len(n)] <- r # possibly recycling r
919                    } else                        } ## else: result = e1  (is "U" diag)
920                      } else {
921                        r <- callGeneric(e1@x, e2)                        r <- callGeneric(e1@x, e2)
922                    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
923                    return(e1)                        e1@x[] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
924                }                }
925                      e1
926                  } else
927                callGeneric(diag2tT.u(e1,e2, "d"), e2)                callGeneric(diag2tT.u(e1,e2, "d"), e2)
928            })            })
929    
930  setMethod("Arith", signature(e1 = "numeric", e2 = "ddiMatrix"),  for(arg1 in c("numeric","logical"))
931    setMethod("Arith", signature(e1 = arg1, e2 = "ddiMatrix"),
932            function(e1,e2) {            function(e1,e2) {
933                n <- e2@Dim[1]; nsq <- n^2                n <- e2@Dim[1]
934                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
935                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
936                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
937                    if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)                    if(e2@diag == "U") {
938                    if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {                        if(any((r <- callGeneric(e1, 1)) != 1)) {
939                        e2@diag <- "N"                        e2@diag <- "N"
940                        if(L1) r <- rep.int(r, n)                            e2@x[seq_len(n)] <- r # possibly recycling r
941                    } else                        } ## else: result = e2  (is "U" diag)
942                      } else {
943                        r <- callGeneric(e1, e2@x)                        r <- callGeneric(e1, e2@x)
944                    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
945                    return(e2)                        e2@x[] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
946                }                }
947                      e2
948                  } else
949                callGeneric(e1, diag2tT.u(e2,e1, "d"))                callGeneric(e1, diag2tT.u(e2,e1, "d"))
950            })            })
951    
952  ## ldi* Arith --> result numeric, potentially ddiMatrix  ## ldi* Arith --> result numeric, potentially ddiMatrix
953  setMethod("Arith", signature(e1 = "ldiMatrix", e2 = "numeric"),  for(arg2 in c("numeric","logical"))
954    setMethod("Arith", signature(e1 = "ldiMatrix", e2 = arg2),
955            function(e1,e2) {            function(e1,e2) {
956                n <- e1@Dim[1]; nsq <- n^2                n <- e1@Dim[1]
957                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
958                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
959                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
960                    if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)                    E <- copyClass(e1, "ddiMatrix", c("diag", "Dim", "Dimnames"))#FIXME: if ok, check=FALSE
961                    if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {                    ## E <- copyClass(e1, "ddiMatrix", check=FALSE)
962                        e1@diag <- "N"                    ## storage.mode(E@x) <- "double"
963                        if(L1) r <- rep.int(r, n)                    if(e1@diag == "U") {
964                    } else                        if(any((r <- callGeneric(1, e2)) != 1)) {
965                              E@diag <- "N"
966                              E@x[seq_len(n)] <- r # possibly recycling r
967                          } ## else: result = E  (is "U" diag)
968                      } else {
969                        r <- callGeneric(e1@x, e2)                        r <- callGeneric(e1@x, e2)
970                    e1 <- copyClass(e1, "ddiMatrix", c("diag", "Dim", "Dimnames"))                        ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
971                    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)  
972                }                }
973                callGeneric(diag2tT.u(e1,e2, "d"), e2)                    E
974                  } else
975                      callGeneric(diag2tT.u(e1,e2, "l"), e2)
976            })            })
977    
978  setMethod("Arith", signature(e1 = "numeric", e2 = "ldiMatrix"),  for(arg1 in c("numeric","logical"))
979    setMethod("Arith", signature(e1 = arg1, e2 = "ldiMatrix"),
980            function(e1,e2) {            function(e1,e2) {
981                n <- e2@Dim[1]; nsq <- n^2                n <- e2@Dim[1]
982                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
983                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
984                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
985                    if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)                    E <- copyClass(e2, "ddiMatrix", c("diag", "Dim", "Dimnames"))#FIXME: if ok, check=FALSE
986                    if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {                    ## E <- copyClass(e2, "ddiMatrix", check=FALSE)
987                        e2@diag <- "N"                    ## storage.mode(E@x) <- "double"
988                        if(L1) r <- rep.int(r, n)                    if(e2@diag == "U") {
989                    } else                        if(any((r <- callGeneric(e1, 1)) != 1)) {
990                              E@diag <- "N"
991                              E@x[seq_len(n)] <- r # possibly recycling r
992                          } ## else: result = E  (is "U" diag)
993                      } else {
994                        r <- callGeneric(e1, e2@x)                        r <- callGeneric(e1, e2@x)
995                    e2 <- copyClass(e2, "ddiMatrix", c("diag", "Dim", "Dimnames"))                        ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
996                    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)  
997                }                }
998                callGeneric(e1, diag2tT.u(e2,e1, "d"))                    E
999                  } else
1000                      callGeneric(e1, diag2tT.u(e2,e1, "l"))
1001            })            })
1002    
1003  ## ddi*: for "Ops" without Arith --> result logical, potentially ldi  ## ddi*: for "Ops" without "Arith": <Compare> or <Logic> --> result logical, potentially ldi
1004  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "numeric"),  ##
1005    ## Note that  ("numeric", "ddiMatrix")  is simply swapped, e.g.,
1006    if(FALSE) {
1007        selectMethod("<", c("numeric","lMatrix"))# Compare
1008        selectMethod("&", c("numeric","lMatrix"))# Logic
1009    } ## so we don't need to define a method here :
1010    
1011    for(arg2 in c("numeric","logical"))
1012    setMethod("Ops", signature(e1 = "ddiMatrix", e2 = arg2),
1013            function(e1,e2) {            function(e1,e2) {
1014                n <- e1@Dim[1]; nsq <- n^2                n <- e1@Dim[1]
1015                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
1016                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
1017                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
1018                    if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)                    E <- copyClass(e1, "ldiMatrix", c("diag", "Dim", "Dimnames"))#FIXME: if ok, check=FALSE
1019                    if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {                    ## E <- copyClass(e1, "ldiMatrix", check=FALSE)
1020                        e1@diag <- "N"                    ## storage.mode(E@x) <- "logical"
1021                        if(L1) r <- rep.int(r, n)                    if(e1@diag == "U") {
1022                    } else                        if(any((r <- callGeneric(1, e2)) != 1)) {
1023                              E@diag <- "N"
1024                              E@x[seq_len(n)] <- r # possibly recycling r
1025                          } ## else: result = E  (is "U" diag)
1026                      } else {
1027                        r <- callGeneric(e1@x, e2)                        r <- callGeneric(e1@x, e2)
1028                    e1 <- copyClass(e1, "ldiMatrix", c("diag", "Dim", "Dimnames"))                        ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
1029                    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)  
1030                }                }
1031                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)  
1032                    } else                    } else
1033                        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"))  
1034            })            })
1035    
1036  ## ldi*: for "Ops" without Arith --> result logical, potentially ldi  ## ldi*: for "Ops" without "Arith": <Compare> or <Logic> --> result logical, potentially ldi
1037  setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "numeric"),  for(arg2 in c("numeric","logical"))
1038    setMethod("Ops", signature(e1 = "ldiMatrix", e2 = arg2),
1039            function(e1,e2) {            function(e1,e2) {
1040                n <- e1@Dim[1]; nsq <- n^2                n <- e1@Dim[1]
1041                f0 <- callGeneric(FALSE, e2)                f0 <- callGeneric(FALSE, e2)
1042                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
1043                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
1044                    if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)  
1045                    if(e1@diag == "U" && any((r <- callGeneric(TRUE, e2)) != 1)) {                    if(e1@diag == "U") {
1046                          if(any((r <- callGeneric(TRUE, e2)) != 1)) {
1047                        e1@diag <- "N"                        e1@diag <- "N"
1048                        if(L1) r <- rep.int(r, n)                            e1@x[seq_len(n)] <- r # possibly recycling r
1049                    } else                        } ## else: result = e1  (is "U" diag)
1050                      } else {
1051                        r <- callGeneric(e1@x, e2)                        r <- callGeneric(e1@x, e2)
1052                    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
1053                    return(e1)                        e1@x[] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
1054                }                }
1055                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)  
1056                    } else                    } else
1057                        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"))  
1058            })            })
1059    
1060    
   
1061  ## Not {"sparseMatrix", "numeric} :  {"denseMatrix", "matrix", ... }  ## Not {"sparseMatrix", "numeric} :  {"denseMatrix", "matrix", ... }
1062  for(other in c("ANY", "Matrix", "dMatrix")) {  for(other in c("ANY", "Matrix", "dMatrix")) {
1063      ## ddi*:      ## ddi*:
1064      setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other),      setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other),
1065                function(e1,e2) callGeneric(diag2tT.u(e1,e2, "d"), e2))                function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2))
1066      setMethod("Ops", signature(e1 = other, e2 = "ddiMatrix"),      setMethod("Ops", signature(e1 = other, e2 = "ddiMatrix"),
1067                function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "d")))                function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "d")))
1068      ## ldi*:      ## ldi*:
1069      setMethod("Ops", signature(e1 = "ldiMatrix", e2 = other),      setMethod("Ops", signature(e1 = "ldiMatrix", e2 = other),
1070                function(e1,e2) callGeneric(diag2tT.u(e1,e2, "l"), e2))                function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2))
1071      setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),      setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),
1072                function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "l")))                function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "l")))
1073  }  }
1074  ## This should *not* dispatch to <dense> methods (in ./Ops.R ), as  
1075  ##  FALSE & <anything> |-> FALSE : hence result should be diagonal:  ## Direct subclasses of "denseMatrix": currently ddenseMatrix, ldense... :
1076  for(cl in diCls) {  dense.subCl <- local({ dM.scl <- getClass("denseMatrix")@subclasses
1077      setMethod("&", signature(e1 = cl, e2 = "ANY"),                         names(dM.scl)[vapply(dM.scl, slot, 0, "distance") == 1] })
1078                function(e1,e2) e1 & as(e2,"Matrix"))  for(DI in diCls) {
1079      setMethod("&", signature(e1 = "ANY", e2 = cl),      dMeth <- if(extends(DI, "dMatrix"))
1080                function(e1,e2) as(e1,"Matrix") & e2)          function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2)
1081      for(c2 in c("denseMatrix", "Matrix")) {      else # "lMatrix", the only other kind for now
1082          setMethod("&", signature(e1 = cl, e2 = c2),          function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2)
1083                    function(e1,e2) e1 & Diagonal(x = diag(e2)))      for(c2 in c(dense.subCl, "Matrix")) {
1084          setMethod("&", signature(e1 = c2, e2 = cl),          for(Fun in c("*", "&")) {
1085                    function(e1,e2) Diagonal(x = diag(e1)) & e2)              setMethod(Fun, signature(e1 = DI, e2 = c2),
1086                          function(e1,e2) callGeneric(e1, Diagonal(x = diag(e2))))
1087                setMethod(Fun, signature(e1 = c2, e2 = DI),
1088                          function(e1,e2) callGeneric(Diagonal(x = diag(e1)), e2))
1089            }
1090            setMethod("^", signature(e1 = c2, e2 = DI),
1091                      function(e1,e2) callGeneric(Diagonal(x = diag(e1)), e2))
1092            for(Fun in c("%%", "%/%", "/")) ## 0 <op> 0 |--> NaN  for these.
1093                setMethod(Fun, signature(e1 = DI, e2 = c2), dMeth)
1094      }      }
1095  }  }
1096    
1097    
1098  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"
1099  ### ----------  any, all: separately here  ### ----------   the last 4: separately here
1100  for(cl in diCls) {  for(cl in diCls) {
1101  setMethod("any", cl,  setMethod("any", cl,
1102            function (x, ..., na.rm) {            function (x, ..., na.rm) {
1103                if(any(x@Dim == 0)) FALSE                if(any(x@Dim == 0)) FALSE
1104                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)
1105            })            })
1106  setMethod("all",  cl, function (x, ..., na.rm) any(x@Dim == 0))  setMethod("all",  cl, function (x, ..., na.rm) {
1107  setMethod("prod", cl, function (x, ..., na.rm) as.numeric(any(x@Dim == 0)))      n <- x@Dim[1]
1108        if(n >= 2) FALSE
1109        else if(n == 0 || x@diag == "U") TRUE
1110        else all(x@x, ..., na.rm = na.rm)
1111    })
1112    setMethod("prod", cl, function (x, ..., na.rm) {
1113        n <- x@Dim[1]
1114        if(n >= 2) 0
1115        else if(n == 0 || x@diag == "U") 1
1116        else ## n == 1, diag = "N" :
1117            prod(x@x, ..., na.rm = na.rm)
1118    })
1119    
1120  setMethod("sum", cl,  setMethod("sum", cl,
1121            function(x, ..., na.rm) {            function(x, ..., na.rm) {
# Line 979  Line 1172 
1172                    invisible(object)                    invisible(object)
1173                }                }
1174            })            })
1175    
1176    rm(dense.subCl, diCls)# not used elsewhere
1177    
1178    setMethod("summary", signature(object = "diagonalMatrix"),
1179              function(object, ...) {
1180                  d <- dim(object)
1181                  r <- summary(object@x, ...)
1182                  attr(r, "header") <-
1183                      sprintf('%d x %d diagonal Matrix of class "%s"',
1184                              d[1], d[2], class(object))
1185                  ## use ole' S3 technology for such a simple case
1186                  class(r) <- c("diagSummary", class(r))
1187                  r
1188              })
1189    
1190    print.diagSummary <- function (x, ...) {
1191        cat(attr(x, "header"),"\n")
1192        class(x) <- class(x)[-1]
1193        print(x, ...)
1194        invisible(x)
1195    }

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

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