SCM

SCM Repository

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

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

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

pkg/R/diagMatrix.R revision 2197, Mon Jun 2 14:34:42 2008 UTC pkg/Matrix/R/diagMatrix.R revision 2863, Mon Dec 31 20:21:11 2012 UTC
# Line 5  Line 5 
5  ##          but *not* diag() extractor!  ##          but *not* diag() extractor!
6  Diagonal <- function(n, x = NULL)  Diagonal <- function(n, x = NULL)
7  {  {
8      ## Allow  Diagonal(4)  and  Diagonal(x=1:5)      ## Allow  Diagonal(4), Diagonal(x=1:5), and  Diagonal(4, TRUE)
9      if(missing(n))      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  ## Pkg 'spdep' had (relatively slow) versions of this as_dsCMatrix_I()  .sparseDiagonal <- function(n, x = 1, uplo = "U",
38  .symDiagonal <- function(n, x = rep.int(1,n), uplo = "U") {                              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      cls <-      m <- length(cols)
46          if(is.double(x)) "dsCMatrix"      if(missing(kind))
47          else if(is.logical(x)) "lsCMatrix"          kind <-
48                if(is.double(x)) "d"
49                else if(is.logical(x)) "l"
50          else { ## for now          else { ## for now
51              storage.mode(x) <- "double"              storage.mode(x) <- "double"
52              "dsCMatrix"                  "d"
53                }
54        else stopifnot(any(kind == c("d","l","n")))
55        stopifnot(is.character(shape), nchar(shape) == 1,
56                  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          }          }
     new(cls, Dim = c(n,n), x = x, uplo = uplo,  
         i = if(n) 0:(n - 1L) else integer(0), p = 0:n)  
79  }  }
80    
81  ### This is modified from a post of Bert Gunter to R-help on  1 Sep 2005.  ## Pkg 'spdep' had (relatively slow) versions of this as_dsCMatrix_I()
82  ### Bert's code built on a post by Andy Liaw who most probably was influenced  .symDiagonal <- function(n, x = rep.int(1,n), uplo = "U")
83  ### by earlier posts, notably one by Scott Chasalow on S-news, 16 Jan 2002      .sparseDiagonal(n, x, uplo, shape = "s")
 ### who posted his bdiag() function written in December 1995.  
84    
85  bdiag <- function(...) {  # instead of   diagU2N(as(Diagonal(n), "CsparseMatrix")), diag = "N" in any case:
86      if(nargs() == 0) return(new("dgCMatrix"))  .trDiagonal <- function(n, x = 1, uplo = "U", unitri=TRUE)
87      ## else :      .sparseDiagonal(n, x, uplo, shape = "t", unitri=unitri)
88      mlist <- if (nargs() == 1) as.list(...) else list(...)  
89      dims <- sapply(mlist, dim)  
90    ## 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
92    ## by earlier posts, notably one by Scott Chasalow on S-news, 16 Jan 2002
93    ## who posted his bdiag() function written in December 1995.
94    if(FALSE)##--- no longer used:
95    .bdiag <- function(lst) {
96        ## block-diagonal matrix [a dgTMatrix] from list of matrices
97        stopifnot(is.list(lst), length(lst) >= 1)
98        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")
102      csdim <- rbind(rep.int(0L, 2),      csdim <- rbind(rep.int(0L, 2),
103                     apply(sapply(mlist, dim), 1, cumsum))                     apply(dims, 1, cumsum))
104      ret <- new("dgTMatrix", Dim = as.integer(csdim[nrow(csdim),]))      r <- new("dgTMatrix")
105        r@Dim <- as.integer(csdim[nrow(csdim),])
106      add1 <- matrix(1:0, 2,2)      add1 <- matrix(1:0, 2,2)
107      for(i in seq_along(mlist)) {      for(i in seq_along(lst)) {
108          indx <- apply(csdim[i:(i+1),] + add1, 2, function(n) n[1]:n[2])          indx <- apply(csdim[i:(i+1),] + add1, 2, function(n) n[1]:n[2])
109          if(is.null(dim(indx))) ## non-square matrix          if(is.null(dim(indx))) ## non-square matrix
110              ret[indx[[1]],indx[[2]]] <- mlist[[i]]              r[indx[[1]],indx[[2]]] <- lst[[i]]
111          else ## square matrix          else ## square matrix
112              ret[indx[,1],indx[,2]] <- mlist[[i]]              r[indx[,1], indx[,2]] <- lst[[i]]
113        }
114        r
115    }
116    ## expand(<mer>) needed something like bdiag() for lower-triangular
117    ## (Tsparse) Matrices; hence Doug Bates provided a much more efficient
118    ##  implementation for those; now extended and generalized:
119    .bdiag <- function(lst) {
120        ## block-diagonal matrix [a dgTMatrix] from list of matrices
121        stopifnot(is.list(lst), (nl <- length(lst)) >= 1)
122    
123        Tlst <- lapply(lapply(lst, Matrix:::as_Csp2), # includes "diagU2N"
124                       as, "TsparseMatrix")
125        if(nl == 1) return(Tlst[[1]])
126        ## else
127        i_off <- c(0L, cumsum(vapply(Tlst, nrow, 1L)))
128        j_off <- c(0L, cumsum(vapply(Tlst, ncol, 1L)))
129    
130        clss <- vapply(Tlst, class, "")
131        typ <- substr(clss, 2, 2)
132        knd <- substr(clss, 1, 1)
133        sym <- typ == "s" # symmetric ones
134        tri <- typ == "t" # triangular ones
135        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")
138            knd [is.n] <- "l"
139        }
140        use.l <- !use.n && all(knd == "l")
141        if(all(sym)) { ## result should be *symmetric*
142            uplos <- vapply(Tlst, slot, ".", "uplo") ## either "U" or "L"
143            tLU <- table(uplos)# of length 1 or 2 ..
144            if(length(tLU) == 1) { ## all "U" or all "L"
145                useU <- uplos[1] == "U"
146            } else { ## length(tLU) == 2, counting "L" and "U"
147                useU <- diff(tLU) >= 0
148                if(useU && (hasL <- tLU[1] > 0))
149                    Tlst[hasL] <- lapply(Tlst[hasL], t)
150                else if(!useU && (hasU <- tLU[2] > 0))
151                    Tlst[hasU] <- lapply(Tlst[hasU], t)
152            }
153            if(use.n) { ## return nsparseMatrix :
154                r <- new("nsTMatrix")
155            } else {
156                r <- new(paste0(if(use.l) "l" else "d", "sTMatrix"))
157                r@x <- unlist(lapply(Tlst, slot, "x"))
158            }
159            r@uplo <- if(useU) "U" else "L"
160        }
161        else if(all(tri) && { ULs <- vapply(Tlst, slot, ".", "uplo")##  "U" or "L"
162                              all(ULs[1L] == ULs[-1L]) } ## all upper or all lower
163           ){ ## *triangular* result
164    
165            if(use.n) { ## return nsparseMatrix :
166                r <- new("ntTMatrix")
167            } else {
168                r <- new(paste0(if(use.l) "l" else "d", "tTMatrix"))
169                r@x <- unlist(lapply(Tlst, slot, "x"))
170            }
171            r@uplo <- ULs[1L]
172        }
173        else {
174            if(any(sym))
175                Tlst[sym] <- lapply(Tlst[sym], as, "generalMatrix")
176            if(use.n) { ## return nsparseMatrix :
177                r <- new("ngTMatrix")
178            } else {
179                r <- new(paste0(if(use.l) "l" else "d", "gTMatrix"))
180                r@x <- unlist(lapply(Tlst, slot, "x"))
181      }      }
182      ## slightly debatable if we really should return Csparse.. :      }
183      as(ret, "CsparseMatrix")      r@Dim <- c(i_off[nl+1], j_off[nl + 1])
184        r@i <- unlist(lapply(1:nl, function(k) Tlst[[k]]@i + i_off[k]))
185        r@j <- unlist(lapply(1:nl, function(k) Tlst[[k]]@j + j_off[k]))
186        r
187    }
188    
189    bdiag <- function(...) {
190        if((nA <- nargs()) == 0) return(new("dgCMatrix"))
191        if(nA == 1 && !is.list(...))
192            return(as(..., "CsparseMatrix"))
193        alis <- if(nA == 1 && is.list(..1)) ..1 else list(...)
194        if(length(alis) == 1)
195            return(as(alis[[1]], "CsparseMatrix"))
196    
197        ## else : two or more arguments
198        as(.bdiag(alis), "CsparseMatrix")
199  }  }
200    
201    
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 91  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 100  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 116  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 124  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 135  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 152  Line 283 
283                Dim = from@Dim, Dimnames = from@Dimnames)                Dim = from@Dim, Dimnames = from@Dimnames)
284        })        })
285    
286    setAs("diagonalMatrix", "nsparseMatrix", function(from) as(from, "nMatrix"))
287    
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 173  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 181  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) 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
320    
321  .diag.2N <- function(m) {  .diag.2N <- function(m) {
# Line 201  Line 336 
336            d <- dim(from)            d <- dim(from)
337            if(d[1] != (n <- d[2])) stop("non-square matrix")            if(d[1] != (n <- d[2])) stop("non-square matrix")
338            if(any(from[row(from) != col(from)] != 0))            if(any(from[row(from) != col(from)] != 0))
339                stop("has non-zero off-diagonal entries")                stop("matrix with non-zero off-diagonals cannot be coerced to \"diagonalMatrix\"")
340            x <- diag(from)            x <- diag(from)
341            if(is.logical(x)) {            if(is.logical(x)) {
342                cl <- "ldiMatrix"                cl <- "ldiMatrix"
343                uni <- all(x)                uni <- allTrue(x) ## uni := {is it unit-diagonal ?}
344            } else {            } else {
345                cl <- "ddiMatrix"                cl <- "ddiMatrix"
346                uni <- all(x == 1)                uni <- allTrue(x == 1)
347                storage.mode(x) <- "double"                storage.mode(x) <- "double"
348            } ## TODO: complex            } ## TODO: complex
349            new(cl, Dim = c(n,n), diag = if(uni) "U" else "N",            new(cl, Dim = c(n,n), diag = if(uni) "U" else "N",
# Line 225  Line 360 
360            x <- diag(from)            x <- diag(from)
361            if(is.logical(x)) {            if(is.logical(x)) {
362                cl <- "ldiMatrix"                cl <- "ldiMatrix"
363                uni <- all(x)                uni <- allTrue(x)
364            } else {            } else {
365                cl <- "ddiMatrix"                cl <- "ddiMatrix"
366                uni <- all(x == 1)                uni <- allTrue(x == 1)
367                storage.mode(x) <- "double"                storage.mode(x) <- "double"
368            }            } ## TODO: complex
369            new(cl, Dim = c(n,n), diag = if(uni) "U" else "N",            new(cl, Dim = c(n,n), diag = if(uni) "U" else "N",
370                x = if(uni) x[FALSE] else x)                x = if(uni) x[FALSE] else x)
371        })        })
# Line 240  Line 375 
375            function(x = 1, nrow, ncol) .diag.x(x))            function(x = 1, nrow, ncol) .diag.x(x))
376    
377  subDiag <- function(x, i, j, ..., drop) {  subDiag <- function(x, i, j, ..., drop) {
378      x <- as(x, "sparseMatrix")      x <- as(x, "CsparseMatrix") ## << was "TsparseMatrix" (Csparse is faster now)
379      x <- if(missing(i))      x <- if(missing(i))
380          x[, j, drop=drop]          x[, j, drop=drop]
381      else if(missing(j))      else if(missing(j))
382          x[i, , drop=drop]          if(nargs() == 4) x[i, , drop=drop] else x[i, drop=drop]
383      else      else
384          x[i,j, drop=drop]          x[i,j, drop=drop]
385      if(isS4(x) && isDiagonal(x)) as(x, "diagonalMatrix") else x      if(isS4(x) && isDiagonal(x)) as(x, "diagonalMatrix") else x
# Line 254  Line 389 
389                           j = "index", drop = "logical"), subDiag)                           j = "index", drop = "logical"), subDiag)
390  setMethod("[", signature(x = "diagonalMatrix", i = "index",  setMethod("[", signature(x = "diagonalMatrix", i = "index",
391                          j = "missing", drop = "logical"),                          j = "missing", drop = "logical"),
392            function(x, i, j, ..., drop) subDiag(x, i=i, drop=drop))            function(x, i, j, ..., drop) {
393                  na <- nargs()
394                  Matrix.msg("diag[i,m,l] : nargs()=", na, .M.level = 2)
395                  if(na == 4)
396                       subDiag(x, i=i, , drop=drop)
397                  else subDiag(x, i=i,   drop=drop)
398              })
399  setMethod("[", signature(x = "diagonalMatrix", i = "missing",  setMethod("[", signature(x = "diagonalMatrix", i = "missing",
400                           j = "index", drop = "logical"),                           j = "index", drop = "logical"),
401            function(x, i, j, ..., drop) subDiag(x, j=j, drop=drop))            function(x, i, j, ..., drop) subDiag(x, j=j, drop=drop))
# Line 264  Line 405 
405  ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch  ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch
406  ## Only(?) current bug:  x[i] <- value  is wrong when  i is *vector*  ## Only(?) current bug:  x[i] <- value  is wrong when  i is *vector*
407  replDiag <- function(x, i, j, ..., value) {  replDiag <- function(x, i, j, ..., value) {
408      x <- as(x, "sparseMatrix")      x <- as(x, "CsparseMatrix")# was "Tsparse.." till 2012-07
409      if(missing(i))      if(missing(i))
410          x[, j] <- value          x[, j] <- value
411      else if(missing(j)) { ##  x[i , ] <- v  *OR*   x[i] <- v      else if(missing(j)) { ##  x[i , ] <- v  *OR*   x[i] <- v
# Line 274  Line 415 
415              x[i, ] <- value              x[i, ] <- value
416          else if(na == 3)          else if(na == 3)
417              x[i] <- value              x[i] <- value
418          else stop("Internal bug: nargs()=",na,"; please report")          else stop(gettextf("Internal bug: nargs()=%d; please report",
419                               na), domain=NA)
420      } else      } else
421          x[i,j] <- value          x[i,j] <- value
422      if(isDiagonal(x)) as(x, "diagonalMatrix") else x      if(isDiagonal(x)) as(x, "diagonalMatrix") else x
# Line 304  Line 446 
446                                   if(any(value != one | is.na(value))) {                                   if(any(value != one | is.na(value))) {
447                                       x@diag <- "N"                                       x@diag <- "N"
448                                       x@x <- rep.int(one, x@Dim[1])                                       x@x <- rep.int(one, x@Dim[1])
449                                   }                                   } else return(x)
450                               }                               }
451                               x@x[ii] <- value                               x@x[ii] <- value
452                               x                               x
453                           } else { ## no longer diagonal, but remain sparse:                           } else { ## no longer diagonal, but remain sparse:
454                               x <- as(x, "sparseMatrix")                               x <- as(x, "TsparseMatrix")
455                               x[i] <- value                               x[i] <- value
456                               x                               x
457                           }                           }
# Line 325  Line 467 
467                                  j = "index", value = "replValue"),                                  j = "index", value = "replValue"),
468                   function(x,i,j, ..., value) replDiag(x, j=j, value=value))                   function(x,i,j, ..., value) replDiag(x, j=j, value=value))
469    
470    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",
471                                    value = "sparseMatrix"),
472                     function (x, i, j, ..., value)
473                     callGeneric(x=x, , j=j, value = as(value, "sparseVector")))
474    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "missing",
475                                    value = "sparseMatrix"),
476                     function (x, i, j, ..., value)
477                     callGeneric(x=x, i=i, , value = as(value, "sparseVector")))
478    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "index",
479                                    value = "sparseMatrix"),
480                     function (x, i, j, ..., value)
481                     callGeneric(x=x, i=i, j=j, value = as(value, "sparseVector")))
482    
483    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",
484                                    value = "sparseVector"),
485                     replDiag)
486    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "missing",
487                                    value = "sparseVector"),
488                     replDiag)
489    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "index",
490                                    value = "sparseVector"),
491                     replDiag)
492    
493    
494  setMethod("t", signature(x = "diagonalMatrix"),  setMethod("t", signature(x = "diagonalMatrix"),
495            function(x) { x@Dimnames <- x@Dimnames[2:1] ; x })            function(x) { x@Dimnames <- x@Dimnames[2:1] ; x })
# Line 371  Line 536 
536  ##       ---------------------  ##       ---------------------
537  ## Note that "ldi" logical are treated as numeric  ## Note that "ldi" logical are treated as numeric
538  diagdiagprod <- function(x, y) {  diagdiagprod <- function(x, y) {
539      n <- dimCheck(x,y)[1]      dimCheck(x,y)
540      if(x@diag != "U") {      if(x@diag != "U") {
541          if(y@diag != "U") {          if(y@diag != "U") {
542              nx <- x@x * y@x              nx <- x@x * y@x
# Line 403  Line 568 
568      dx <- dim(x)      dx <- dim(x)
569      dy <- dim(y)      dy <- dim(y)
570      if(dx[2] != dy[1]) stop("non-matching dimensions")      if(dx[2] != dy[1]) stop("non-matching dimensions")
     n <- dx[1]  
571      as(if(x@diag == "U") y else x@x * y, "Matrix")      as(if(x@diag == "U") y else x@x * y, "Matrix")
572  }  }
573  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),
# Line 472  Line 636 
636  ##        function(x, y = NULL) {  ##        function(x, y = NULL) {
637  ##           })  ##           })
638    
639  ## setMethod("tcrossprod", signature(x = "denseMatrix", y = "diagonalMatrix"),  Cspdiagprod <- function(x, y) {
640  ##        function(x, y = NULL) {      dx <- dim(x <- .Call(Csparse_diagU2N, x))
641  ##           })      dy <- dim(y)
642        if(dx[2] != dy[1]) stop("non-matching dimensions")
643        if(y@diag == "N") {
644            if(!all(y@x[1L] == y@x[-1L]) && is(x, "symmetricMatrix"))
645                x <- as(x, "generalMatrix")
646            ind <- rep.int(seq_len(dx[2]), x@p[-1] - x@p[-dx[2]-1L])
647            x@x <- x@x * y@x[ind]
648        }
649        if(is(x, "compMatrix") && length(xf <- x@factors)) {
650            ## instead of dropping all factors, be smart about some
651            ## TODO ......
652            x@factors <- list()
653        }
654        x
655    }
656    
657    diagCspprod <- function(x, y) {
658        dx <- dim(x)
659        dy <- dim(y <- .Call(Csparse_diagU2N, y))
660        if(dx[2] != dy[1]) stop("non-matching dimensions")
661        if(x@diag == "N") {
662            if(!all(x@x[1L] == x@x[-1L]) && is(y, "symmetricMatrix"))
663                y <- as(y, "generalMatrix")
664            y@x <- y@x * x@x[y@i + 1L]
665        }
666        if(is(y, "compMatrix") && length(yf <- y@factors)) {
667            ## instead of dropping all factors, be smart about some
668            ## TODO
669            keep <- character()
670            if(iLU <- names(yf) == "LU") {
671                ## TODO keep <- "LU"
672            }
673            y@factors <- yf[keep]
674        }
675        y
676    }
677    
678    setMethod("crossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
679              function(x, y = NULL) diagCspprod(x, y))
680    
681  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
682            function(x, y = NULL) crossprod(as(x, "sparseMatrix"), y))            function(x, y = NULL) diagCspprod(x, as(y, "CsparseMatrix")))
683    
684    ## Prefer calling diagCspprod to Cspdiagprod if going to transpose anyway
685    ##  x'y == (y'x)'
686    setMethod("crossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
687              function(x, y = NULL) t(diagCspprod(y, x)))
688    
689  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
690            function(x, y = NULL) crossprod(x, as(y, "sparseMatrix")))            function(x, y = NULL) t(diagCspprod(y, as(x, "Csparsematrix"))))
691    
692    setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
693              function(x, y = NULL) diagCspprod(x, t(y)))
694    
695  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
696            function(x, y = NULL) tcrossprod(as(x, "sparseMatrix"), y))            function(x, y = NULL) diagCspprod(x, t(as(y, "CsparseMatrix"))))
697    
698    setMethod("tcrossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
699              function(x, y = NULL) Cspdiagprod(x, y))
700    
701  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
702            function(x, y = NULL) tcrossprod(x, as(y, "sparseMatrix")))            function(x, y = NULL) Cspdiagprod(as(x, "CsparseMatrix"), y))
703    
704    setMethod("%*%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
705              function(x, y) diagCspprod(x, y))
706    
 ## FIXME?: In theory, this can be done *FASTER*, in some cases, via tapply1()  
707  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
708            function(x, y) as(x, "sparseMatrix") %*% y)            function(x, y) diagCspprod(as(x, "CsparseMatrix"), y))
709    
710  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
711            function(x, y) x %*% as(y, "sparseMatrix"))            function(x, y) Cspdiagprod(as(x, "CsparseMatrix"), y))
712  ## 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)  
713  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
714            function(x, y) x %*% as(y, "CsparseMatrix"))            function(x, y) Cspdiagprod(x, y))
715  ## NB: this is *not* needed for Tsparse & Rsparse  
716  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*
717  ##       do indeed work by going through sparse (and *not* ddense)!  ##       do indeed work by going through sparse (and *not* ddense)!
718    
   
   
719  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),
720            function(a, b, ...) {            function(a, b, ...) {
721                a@x <- 1/ a@x                a@x <- 1/ a@x
# Line 515  Line 724 
724            })            })
725    
726  solveDiag <- function(a, b, ...) {  solveDiag <- function(a, b, ...) {
727      if((n <- a@Dim[1]) != nrow(b))      if(a@Dim[1] != nrow(b))
728          stop("incompatible matrix dimensions")          stop("incompatible matrix dimensions")
729      ## trivially invert a 'in place' and multiply:      ## trivially invert a 'in place' and multiply:
730      a@x <- 1/ a@x      a@x <- 1/ a@x
# Line 534  Line 743 
743  ###---------------- <Ops> (<Arith>, <Logic>, <Compare> ) ----------------------  ###---------------- <Ops> (<Arith>, <Logic>, <Compare> ) ----------------------
744    
745  ## 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.)  
746  diagOdiag <- function(e1,e2) {  diagOdiag <- function(e1,e2) {
747      ## result should also be diagonal _ if possible _      ## result should also be diagonal _ if possible _
748      r <- callGeneric(.diag.x(e1), .diag.x(e2)) # error if not "compatible"      r <- callGeneric(.diag.x(e1), .diag.x(e2)) # error if not "compatible"
# Line 542  Line 750 
750      r00 <- callGeneric(if(is.numeric(e1@x)) 0 else FALSE,      r00 <- callGeneric(if(is.numeric(e1@x)) 0 else FALSE,
751                         if(is.numeric(e2@x)) 0 else FALSE)                         if(is.numeric(e2@x)) 0 else FALSE)
752      if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* diagonal      if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* diagonal
753          if(is.numeric(r)) {          if(is.numeric(r)) { # "double" *or* "integer"
754              if(is.numeric(e2@x)) {              if(is.numeric(e2@x)) {
755                  e2@x <- r; return(.diag.2N(e2)) }                  e2@x <- r; return(.diag.2N(e2)) }
756              if(!is.numeric(e1@x))              if(!is.numeric(e1@x))
757                  ## e.g. e1, e2 are logical;                  ## e.g. e1, e2 are logical;
758                  e1 <- as(e1, "dMatrix")                  e1 <- as(e1, "dMatrix")
759                if(!is.double(r)) r <- as.double(r)
760          }          }
761          else if(is.logical(r))          else if(is.logical(r))
762              e1 <- as(e1, "lMatrix")              e1 <- as(e1, "lMatrix")
763          else stop("intermediate 'r' is of type", typeof(r))          else stop(gettextf("intermediate 'r' is of type %s",
764                               typeof(r)), domain=NA)
765          e1@x <- r          e1@x <- r
766          .diag.2N(e1)          .diag.2N(e1)
767      }      }
768      else { ## result not diagonal, but at least symmetric:      else { ## result not diagonal, but at least symmetric:
769            ## e.g., m == m
770          isNum <- (is.numeric(r) || is.numeric(r00))          isNum <- (is.numeric(r) || is.numeric(r00))
771          isLog <- (is.logical(r) || is.logical(r00))          isLog <- (is.logical(r) || is.logical(r00))
772            Matrix.msg("exploding <diag> o <diag> into dense matrix", .M.level = 2)
         if(getOption("verbose"))  
             message("exploding  <diag>  o  <diag>  into dense matrix")  
773          d <- e1@Dim          d <- e1@Dim
774          n <- d[1]          n <- d[1]
775          stopifnot(length(r) == n)          stopifnot(length(r) == n)
776            if(isNum && !is.double(r)) r <- as.double(r)
777          xx <- as.vector(matrix(rbind(r, matrix(r00,n,n)), n,n))          xx <- as.vector(matrix(rbind(r, matrix(r00,n,n)), n,n))
778          newcl <-          newcl <-
779              paste(if(isNum) "d" else if(isLog) {              paste0(if(isNum) "d" else if(isLog) {
780                  if(!any(is.na(r)) && !any(is.na(r00))) "n" else "l"                  if(!any(is.na(r)) && !any(is.na(r00))) "n" else "l"
781              } else stop("not yet implemented .. please report")              } else stop("not yet implemented .. please report"), "syMatrix")
                   ,  
                   "syMatrix", sep='')  
782    
783          new(newcl, Dim = e1@Dim, Dimnames = e1@Dimnames, x = xx)          new(newcl, Dim = e1@Dim, Dimnames = e1@Dimnames, x = xx)
784      }      }
# Line 588  Line 796 
796              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)
797  }  }
798    
799  ## FIXME:    diagonal  o  triangular  |-->  triangular  ## diagonal  o  triangular  |-->  triangular
800  ## -----     diagonal  o  symmetric   |-->  symmetric  ## diagonal  o  symmetric   |-->  symmetric
801  ##    {also when other is sparse: do these "here" --  ##    {also when other is sparse: do these "here" --
802  ##     before conversion to sparse, since that loses "diagonality"}  ##     before conversion to sparse, since that loses "diagonality"}
803    diagOtri <- function(e1,e2) {
804        ## result must be triangular
805        r <- callGeneric(d1 <- .diag.x(e1), diag(e2)) # error if not "compatible"
806        ## Check what happens with non-diagonals, i.e. (0 o 0), (FALSE o 0), ...:
807        e1.0 <- if(.n1 <- is.numeric(d1   )) 0 else FALSE
808        r00 <- callGeneric(e1.0, if(.n2 <- is.numeric(e2[0])) 0 else FALSE)
809        if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* triangular
810            diag(e2) <- r
811            ## check what happens "in the triangle"
812            e2.2 <- if(.n2) 2 else TRUE
813            if(!callGeneric(e1.0, e2.2) == e2.2) { # values "in triangle" can change:
814                n <- dim(e2)[1L]
815                it <- indTri(n, upper = (e2@uplo == "U"))
816                e2[it] <- callGeneric(e1.0, e2[it])
817            }
818            e2
819        }
820        else { ## result not triangular ---> general
821            rr <- as(e2, "generalMatrix")
822            diag(rr) <- r
823            rr
824        }
825    }
826    
827    
828    setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "triangularMatrix"),
829              diagOtri)
830    ## For the reverse,  Ops == "Arith" | "Compare" | "Logic"
831    ##   'Arith'  :=  '"+"', '"-"', '"*"', '"^"', '"%%"', '"%/%"', '"/"'
832    setMethod("Arith", signature(e1 = "triangularMatrix", e2 = "diagonalMatrix"),
833              function(e1,e2)
834          { ## this must only trigger for *dense* e1
835              switch(.Generic,
836                     "+" = .Call(dtrMatrix_addDiag, as(e1,"dtrMatrix"),   .diag.x(e2)),
837                     "-" = .Call(dtrMatrix_addDiag, as(e1,"dtrMatrix"), - .diag.x(e2)),
838                     "*" = {
839                         n <- e2@Dim[1L]
840                         d2 <- if(e2@diag == "U") { # unit-diagonal
841                             d <- rep.int(as1(e2@x), n)
842                             e2@x <- d
843                             e2@diag <- "N"
844                             d
845                         } else e2@x
846                         e2@x <- diag(e1) * d2
847                         e2
848                     },
849                     "^" = { ## will be dense ( as  <ANY> ^ 0 == 1 ):
850                         e1 ^ as(e2, "denseMatrix")
851                     },
852                     ## otherwise:
853                     callGeneric(e1, diag2Tsmart(e2,e1)))
854    })
855    
856    ## Compare --> 'swap' (e.g.   e1 < e2   <==>  e2 > e1 ):
857    setMethod("Compare", signature(e1 = "triangularMatrix", e2 = "diagonalMatrix"),
858              .Cmp.swap)
859    ## '&' and "|'  are commutative:
860    setMethod("Logic", signature(e1 = "triangularMatrix", e2 = "diagonalMatrix"),
861              function(e1,e2) callGeneric(e2, e1))
862    
863  ## For almost everything else, diag* shall be treated "as sparse" :  ## For almost everything else, diag* shall be treated "as sparse" :
864  ## These are cheap implementations via coercion  ## These are cheap implementations via coercion
# Line 615  Line 882 
882  ##       Compare --> logical  ##       Compare --> logical
883  ##       Logic   --> logical: "lMatrix"  ##       Logic   --> logical: "lMatrix"
884    
885  ##  other = "numeric" : stay diagonal if possible  ## Other = "numeric" : stay diagonal if possible
886  ## ddi*: Arith: result numeric, potentially ddiMatrix  ## ddi*: Arith: result numeric, potentially ddiMatrix
887  setMethod("Arith", signature(e1 = "ddiMatrix", e2 = "numeric"),  for(arg2 in c("numeric","logical"))
888    setMethod("Arith", signature(e1 = "ddiMatrix", e2 = arg2),
889            function(e1,e2) {            function(e1,e2) {
890                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]
891                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
892                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
893                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
894                    if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)                    if(e1@diag == "U") {
895                    if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {                        if(any((r <- callGeneric(1, e2)) != 1)) {
896                        e1@diag <- "N"                        e1@diag <- "N"
897                        if(L1) r <- rep.int(r, n)                            e1@x[seq_len(n)] <- r # possibly recycling r
898                    } else                        } ## else: result = e1  (is "U" diag)
899                      } else {
900                        r <- callGeneric(e1@x, e2)                        r <- callGeneric(e1@x, e2)
901                    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
902                    return(e1)                        e1@x[] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
903                }                }
904                      e1
905                  } else
906                callGeneric(diag2tT.u(e1,e2, "d"), e2)                callGeneric(diag2tT.u(e1,e2, "d"), e2)
907            })            })
908    
909  setMethod("Arith", signature(e1 = "numeric", e2 = "ddiMatrix"),  for(arg1 in c("numeric","logical"))
910    setMethod("Arith", signature(e1 = arg1, e2 = "ddiMatrix"),
911            function(e1,e2) {            function(e1,e2) {
912                n <- e2@Dim[1]; nsq <- n*n                n <- e2@Dim[1]
913                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
914                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
915                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
916                    if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)                    if(e2@diag == "U") {
917                    if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {                        if(any((r <- callGeneric(e1, 1)) != 1)) {
918                        e2@diag <- "N"                        e2@diag <- "N"
919                        if(L1) r <- rep.int(r, n)                            e2@x[seq_len(n)] <- r # possibly recycling r
920                    } else                        } ## else: result = e2  (is "U" diag)
921                      } else {
922                        r <- callGeneric(e1, e2@x)                        r <- callGeneric(e1, e2@x)
923                    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
924                    return(e2)                        e2@x[] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
925                }                }
926                      e2
927                  } else
928                callGeneric(e1, diag2tT.u(e2,e1, "d"))                callGeneric(e1, diag2tT.u(e2,e1, "d"))
929            })            })
930    
931  ## ldi* Arith --> result numeric, potentially ddiMatrix  ## ldi* Arith --> result numeric, potentially ddiMatrix
932  setMethod("Arith", signature(e1 = "ldiMatrix", e2 = "numeric"),  for(arg2 in c("numeric","logical"))
933    setMethod("Arith", signature(e1 = "ldiMatrix", e2 = arg2),
934            function(e1,e2) {            function(e1,e2) {
935                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]
936                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
937                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
938                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
939                    if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)                    E <- copyClass(e1, "ddiMatrix", c("diag", "Dim", "Dimnames"))#FIXME: if ok, check=FALSE
940                    if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {                    ## E <- copyClass(e1, "ddiMatrix", check=FALSE)
941                        e1@diag <- "N"                    ## storage.mode(E@x) <- "double"
942                        if(L1) r <- rep.int(r, n)                    if(e1@diag == "U") {
943                    } else                        if(any((r <- callGeneric(1, e2)) != 1)) {
944                              E@diag <- "N"
945                              E@x[seq_len(n)] <- r # possibly recycling r
946                          } ## else: result = E  (is "U" diag)
947                      } else {
948                        r <- callGeneric(e1@x, e2)                        r <- callGeneric(e1@x, e2)
949                    e1 <- copyClass(e1, "ddiMatrix", c("diag", "Dim", "Dimnames"))                        ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
950                    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)  
951                }                }
952                callGeneric(diag2tT.u(e1,e2, "d"), e2)                    E
953                  } else
954                      callGeneric(diag2tT.u(e1,e2, "l"), e2)
955            })            })
956    
957  setMethod("Arith", signature(e1 = "numeric", e2 = "ldiMatrix"),  for(arg1 in c("numeric","logical"))
958    setMethod("Arith", signature(e1 = arg1, e2 = "ldiMatrix"),
959            function(e1,e2) {            function(e1,e2) {
960                n <- e2@Dim[1]; nsq <- n*n                n <- e2@Dim[1]
961                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
962                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
963                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
964                    if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)                    E <- copyClass(e2, "ddiMatrix", c("diag", "Dim", "Dimnames"))#FIXME: if ok, check=FALSE
965                    if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {                    ## E <- copyClass(e2, "ddiMatrix", check=FALSE)
966                        e2@diag <- "N"                    ## storage.mode(E@x) <- "double"
967                        if(L1) r <- rep.int(r, n)                    if(e2@diag == "U") {
968                    } else                        if(any((r <- callGeneric(e1, 1)) != 1)) {
969                              E@diag <- "N"
970                              E@x[seq_len(n)] <- r # possibly recycling r
971                          } ## else: result = E  (is "U" diag)
972                      } else {
973                        r <- callGeneric(e1, e2@x)                        r <- callGeneric(e1, e2@x)
974                    e2 <- copyClass(e2, "ddiMatrix", c("diag", "Dim", "Dimnames"))                        ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
975                    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)  
976                }                }
977                callGeneric(e1, diag2tT.u(e2,e1, "d"))                    E
978                  } else
979                      callGeneric(e1, diag2tT.u(e2,e1, "l"))
980            })            })
981    
982  ## ddi*: for "Ops" without Arith --> result logical, potentially ldi  ## ddi*: for "Ops" without "Arith": <Compare> or <Logic> --> result logical, potentially ldi
983  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "numeric"),  ##
984    ## Note that  ("numeric", "ddiMatrix")  is simply swapped, e.g.,
985    if(FALSE) {
986        selectMethod("<", c("numeric","lMatrix"))# Compare
987        selectMethod("&", c("numeric","lMatrix"))# Logic
988    } ## so we don't need to define a method here :
989    
990    for(arg2 in c("numeric","logical"))
991    setMethod("Ops", signature(e1 = "ddiMatrix", e2 = arg2),
992            function(e1,e2) {            function(e1,e2) {
993                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]
994                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
995                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
996                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
997                    if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)                    E <- copyClass(e1, "ldiMatrix", c("diag", "Dim", "Dimnames"))#FIXME: if ok, check=FALSE
998                    if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {                    ## E <- copyClass(e1, "ldiMatrix", check=FALSE)
999                        e1@diag <- "N"                    ## storage.mode(E@x) <- "logical"
1000                        if(L1) r <- rep.int(r, n)                    if(e1@diag == "U") {
1001                    } else                        if(any((r <- callGeneric(1, e2)) != 1)) {
1002                              E@diag <- "N"
1003                              E@x[seq_len(n)] <- r # possibly recycling r
1004                          } ## else: result = E  (is "U" diag)
1005                      } else {
1006                        r <- callGeneric(e1@x, e2)                        r <- callGeneric(e1@x, e2)
1007                    e1 <- copyClass(e1, "ldiMatrix", c("diag", "Dim", "Dimnames"))                        ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
1008                    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)  
1009                }                }
1010                callGeneric(diag2tT.u(e1,e2, "l"), e2)                    E
           })  
   
 setMethod("Ops", signature(e1 = "numeric", e2 = "ddiMatrix"),  
           function(e1,e2) {  
               n <- e2@Dim[1]; nsq <- n*n  
               f0 <- callGeneric(e1, 0)  
               if(all(is0(f0))) { # remain diagonal  
                   L1 <- (le <- length(e1)) == 1L  
                   if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)  
                   if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {  
                       e2@diag <- "N"  
                       if(L1) r <- rep.int(r, n)  
1011                    } else                    } else
1012                        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"))  
1013            })            })
1014    
1015  ## ldi*: for "Ops" without Arith --> result logical, potentially ldi  ## ldi*: for "Ops" without "Arith": <Compare> or <Logic> --> result logical, potentially ldi
1016  setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "numeric"),  for(arg2 in c("numeric","logical"))
1017    setMethod("Ops", signature(e1 = "ldiMatrix", e2 = arg2),
1018            function(e1,e2) {            function(e1,e2) {
1019                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]
1020                f0 <- callGeneric(FALSE, e2)                f0 <- callGeneric(FALSE, e2)
1021                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
1022                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
1023                    if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)  
1024                    if(e1@diag == "U" && any((r <- callGeneric(TRUE, e2)) != 1)) {                    if(e1@diag == "U") {
1025                          if(any((r <- callGeneric(TRUE, e2)) != 1)) {
1026                        e1@diag <- "N"                        e1@diag <- "N"
1027                        if(L1) r <- rep.int(r, n)                            e1@x[seq_len(n)] <- r # possibly recycling r
1028                    } else                        } ## else: result = e1  (is "U" diag)
1029                      } else {
1030                        r <- callGeneric(e1@x, e2)                        r <- callGeneric(e1@x, e2)
1031                    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
1032                    return(e1)                        e1@x[] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
1033                }                }
1034                callGeneric(diag2tT.u(e1,e2, "l"), e2)                    e1
           })  
   
 setMethod("Ops", signature(e1 = "numeric", e2 = "ldiMatrix"),  
           function(e1,e2) {  
               n <- e2@Dim[1]; nsq <- n*n  
               f0 <- callGeneric(e1, FALSE)  
               if(all(is0(f0))) { # remain diagonal  
                   L1 <- (le <- length(e1)) == 1L  
                   if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)  
                   if(e2@diag == "U" && any((r <- callGeneric(e1, TRUE)) != 1)) {  
                       e2@diag <- "N"  
                       if(L1) r <- rep.int(r, n)  
1035                    } else                    } else
1036                        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"))  
1037            })            })
1038    
1039    
   
1040  ## Not {"sparseMatrix", "numeric} :  {"denseMatrix", "matrix", ... }  ## Not {"sparseMatrix", "numeric} :  {"denseMatrix", "matrix", ... }
1041  for(other in c("ANY", "Matrix", "dMatrix")) {  for(other in c("ANY", "Matrix", "dMatrix")) {
1042      ## ddi*:      ## ddi*:
1043      setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other),      setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other),
1044                function(e1,e2) callGeneric(diag2tT.u(e1,e2, "d"), e2))                function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2))
1045      setMethod("Ops", signature(e1 = other, e2 = "ddiMatrix"),      setMethod("Ops", signature(e1 = other, e2 = "ddiMatrix"),
1046                function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "d")))                function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "d")))
1047      ## ldi*:      ## ldi*:
1048      setMethod("Ops", signature(e1 = "ldiMatrix", e2 = other),      setMethod("Ops", signature(e1 = "ldiMatrix", e2 = other),
1049                function(e1,e2) callGeneric(diag2tT.u(e1,e2, "l"), e2))                function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2))
1050      setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),      setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),
1051                function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "l")))                function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "l")))
1052  }  }
1053  ## This should *not* dispatch to <dense> methods (in ./Ops.R ), as  
1054  ##  FALSE & <anything> |-> FALSE : hence result should be diagonal:  ## Direct subclasses of "denseMatrix": currently ddenseMatrix, ldense... :
1055  for(cl in diCls) {  dense.subCl <- local({ dM.scl <- getClass("denseMatrix")@subclasses
1056      setMethod("&", signature(e1 = cl, e2 = "ANY"),                         names(dM.scl)[vapply(dM.scl, slot, 0, "distance") == 1] })
1057                function(e1,e2) e1 & as(e2,"Matrix"))  for(DI in diCls) {
1058      setMethod("&", signature(e1 = "ANY", e2 = cl),      dMeth <- if(extends(DI, "dMatrix"))
1059                function(e1,e2) as(e1,"Matrix") & e2)          function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2)
1060      for(c2 in c("denseMatrix", "Matrix")) {      else # "lMatrix", the only other kind for now
1061          setMethod("&", signature(e1 = cl, e2 = c2),          function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2)
1062                    function(e1,e2) e1 & Diagonal(x = diag(e2)))      for(c2 in c(dense.subCl, "Matrix")) {
1063          setMethod("&", signature(e1 = c2, e2 = cl),          for(Fun in c("*", "&")) {
1064                    function(e1,e2) Diagonal(x = diag(e1)) & e2)              setMethod(Fun, signature(e1 = DI, e2 = c2),
1065                          function(e1,e2) callGeneric(e1, Diagonal(x = diag(e2))))
1066                setMethod(Fun, signature(e1 = c2, e2 = DI),
1067                          function(e1,e2) callGeneric(Diagonal(x = diag(e1)), e2))
1068            }
1069            setMethod("^", signature(e1 = c2, e2 = DI),
1070                      function(e1,e2) callGeneric(Diagonal(x = diag(e1)), e2))
1071            for(Fun in c("%%", "%/%", "/")) ## 0 <op> 0 |--> NaN  for these.
1072                setMethod(Fun, signature(e1 = DI, e2 = c2), dMeth)
1073      }      }
1074  }  }
1075    
1076    
1077  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"
1078  ### ----------  any, all: separately here  ### ----------   the last 4: separately here
1079  for(cl in diCls) {  for(cl in diCls) {
1080  setMethod("any", cl,  setMethod("any", cl,
1081            function (x, ..., na.rm) {            function (x, ..., na.rm) {
1082                if(any(x@Dim == 0)) FALSE                if(any(x@Dim == 0)) FALSE
1083                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)
1084            })            })
1085  setMethod("all",  cl, function (x, ..., na.rm) any(x@Dim == 0))  setMethod("all",  cl, function (x, ..., na.rm) {
1086  setMethod("prod", cl, function (x, ..., na.rm) as.numeric(any(x@Dim == 0)))      n <- x@Dim[1]
1087        if(n >= 2) FALSE
1088        else if(n == 0 || x@diag == "U") TRUE
1089        else all(x@x, ..., na.rm = na.rm)
1090    })
1091    setMethod("prod", cl, function (x, ..., na.rm) {
1092        n <- x@Dim[1]
1093        if(n >= 2) 0
1094        else if(n == 0 || x@diag == "U") 1
1095        else ## n == 1, diag = "N" :
1096            prod(x@x, ..., na.rm = na.rm)
1097    })
1098    
1099  setMethod("sum", cl,  setMethod("sum", cl,
1100            function(x, ..., na.rm) {            function(x, ..., na.rm) {
# Line 847  Line 1133 
1133      invisible(x)      invisible(x)
1134  }  }
1135    
1136    ## somewhat consistent with "print" for sparseMatrix :
1137    setMethod("print", signature(x = "diagonalMatrix"), prDiag)
1138    
1139  setMethod("show", signature(object = "diagonalMatrix"),  setMethod("show", signature(object = "diagonalMatrix"),
1140            function(object) {            function(object) {
1141                d <- dim(object)                d <- dim(object)
1142                cl <- class(object)                cl <- class(object)
1143                cat(sprintf('%d x %d diagonal matrix of class "%s"\n',                cat(sprintf('%d x %d diagonal matrix of class "%s"',
1144                            d[1], d[2], cl))                            d[1], d[2], cl))
1145                  if(d[1] < 50) {
1146                      cat("\n")
1147                prDiag(object)                prDiag(object)
1148                  } else {
1149                      cat(", with diagonal entries\n")
1150                      show(diag(object))
1151                      invisible(object)
1152                  }
1153            })            })
1154    
1155    rm(dense.subCl, diCls)# not used elsewhere
1156    
1157    setMethod("summary", signature(object = "diagonalMatrix"),
1158              function(object, ...) {
1159                  d <- dim(object)
1160                  r <- summary(object@x, ...)
1161                  attr(r, "header") <-
1162                      sprintf('%d x %d diagonal Matrix of class "%s"',
1163                              d[1], d[2], class(object))
1164                  ## use ole' S3 technology for such a simple case
1165                  class(r) <- c("diagSummary", class(r))
1166                  r
1167              })
1168    
1169    print.diagSummary <- function (x, ...) {
1170        cat(attr(x, "header"),"\n")
1171        class(x) <- class(x)[-1]
1172        print(x, ...)
1173        invisible(x)
1174    }

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

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