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 2185, Sat Apr 26 20:33:16 2008 UTC pkg/Matrix/R/diagMatrix.R revision 2904, Tue Sep 10 19:43:53 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  ## 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      ## slightly debatable if we really should return Csparse.. :      r
115      as(ret, "CsparseMatrix")  }
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, 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        }
183        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)                   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 201  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 225  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 240  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 254  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 264  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 274  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 293  Line 436 
436                           replDiag(x, i=i, , value=value)                           replDiag(x, i=i, , value=value)
437                   })                   })
438    
439  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "matrix", # 2-col.matrix  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",
461                                    i = "matrix", # 2-col.matrix
462                                  j = "missing", value = "replValue"),                                  j = "missing", value = "replValue"),
463                   function(x,i,j, ..., value) {                   function(x,i,j, ..., value) {
464                       if(ncol(i) == 2) {                       if(ncol(i) == 2) {
465                           if(all((ii <- i[,1]) == i[,2])) { # replace in diagonal only                           if(all((ii <- i[,1]) == i[,2])) { # replace in diagonal only
466                                 if(x@diag == "U") {
467                                     one <- as1(x@x)
468                                     if(any(value != one | is.na(value))) {
469                                         x@diag <- "N"
470                                         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 313  Line 485 
485                       }                       }
486                   })                   })
487    
488  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing",  
489                                  j = "index", value = "replValue"),  ## value = "sparseMatrix":
490                   function(x,i,j, ..., value) replDiag(x, j=j, value=value))  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",
491                                    value = "sparseMatrix"),
492                     function (x, i, j, ..., value)
493                     callGeneric(x=x, , j=j, value = as(value, "sparseVector")))
494    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "missing",
495                                    value = "sparseMatrix"),
496                     function (x, i, j, ..., value)
497                     callGeneric(x=x, i=i, , value = as(value, "sparseVector")))
498    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "index",
499                                    value = "sparseMatrix"),
500                     function (x, i, j, ..., value)
501                     callGeneric(x=x, i=i, j=j, value = as(value, "sparseVector")))
502    
503    ## value = "sparseVector":
504    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",
505                                    value = "sparseVector"),
506                     replDiag)
507    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "missing",
508                                    value = "sparseVector"),
509                     replDiag)
510    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "index",
511                                    value = "sparseVector"),
512                     replDiag)
513    
514    
515  setMethod("t", signature(x = "diagonalMatrix"),  setMethod("t", signature(x = "diagonalMatrix"),
516            function(x) { x@Dimnames <- x@Dimnames[2:1] ; x })            function(x) { x@Dimnames <- x@Dimnames[2:1] ; x })
517    
518  setMethod("isDiagonal", signature(object = "diagonalMatrix"),  setMethod("isDiagonal",   "diagonalMatrix", function(object) TRUE)
519            function(object) TRUE)  setMethod("isTriangular", "diagonalMatrix", function(object, ...) TRUE)
520  setMethod("isTriangular", signature(object = "diagonalMatrix"),  setMethod("isSymmetric",  "diagonalMatrix", function(object, ...) TRUE)
           function(object) TRUE)  
 setMethod("isSymmetric", signature(object = "diagonalMatrix"),  
           function(object, ...) TRUE)  
521    
522  setMethod("symmpart", signature(x = "diagonalMatrix"), function(x) x)  setMethod("symmpart", signature(x = "diagonalMatrix"), function(x) x)
523  setMethod("skewpart", signature(x = "diagonalMatrix"), setZero)  setMethod("skewpart", signature(x = "diagonalMatrix"), setZero)
# Line 363  Line 554 
554  ##       ---------------------  ##       ---------------------
555  ## Note that "ldi" logical are treated as numeric  ## Note that "ldi" logical are treated as numeric
556  diagdiagprod <- function(x, y) {  diagdiagprod <- function(x, y) {
557      n <- dimCheck(x,y)[1]      dimCheck(x,y)
558      if(x@diag != "U") {      if(x@diag != "U") {
559          if(y@diag != "U") {          if(y@diag != "U") {
560              nx <- x@x * y@x              nx <- x@x * y@x
# Line 395  Line 586 
586      dx <- dim(x)      dx <- dim(x)
587      dy <- dim(y)      dy <- dim(y)
588      if(dx[2] != dy[1]) stop("non-matching dimensions")      if(dx[2] != dy[1]) stop("non-matching dimensions")
     n <- dx[1]  
589      as(if(x@diag == "U") y else x@x * y, "Matrix")      as(if(x@diag == "U") y else x@x * y, "Matrix")
590  }  }
591  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),
# Line 464  Line 654 
654  ##        function(x, y = NULL) {  ##        function(x, y = NULL) {
655  ##           })  ##           })
656    
657  ## setMethod("tcrossprod", signature(x = "denseMatrix", y = "diagonalMatrix"),  Cspdiagprod <- function(x, y) {
658  ##        function(x, y = NULL) {      dx <- dim(x <- .Call(Csparse_diagU2N, x))
659  ##           })      dy <- dim(y)
660        if(dx[2] != dy[1]) stop("non-matching dimensions")
661        if(y@diag == "N") {
662            if(!all(y@x[1L] == y@x[-1L]) && is(x, "symmetricMatrix"))
663                x <- as(x, "generalMatrix")
664            ind <- rep.int(seq_len(dx[2]), x@p[-1] - x@p[-dx[2]-1L])
665            x@x <- x@x * y@x[ind]
666        }
667        if(is(x, "compMatrix") && length(xf <- x@factors)) {
668            ## instead of dropping all factors, be smart about some
669            ## TODO ......
670            x@factors <- list()
671        }
672        x
673    }
674    
675    diagCspprod <- function(x, y) {
676        dx <- dim(x)
677        dy <- dim(y <- .Call(Csparse_diagU2N, y))
678        if(dx[2] != dy[1]) stop("non-matching dimensions")
679        if(x@diag == "N") {
680            if(!all(x@x[1L] == x@x[-1L]) && is(y, "symmetricMatrix"))
681                y <- as(y, "generalMatrix")
682            y@x <- y@x * x@x[y@i + 1L]
683        }
684        if(is(y, "compMatrix") && length(yf <- y@factors)) {
685            ## instead of dropping all factors, be smart about some
686            ## TODO
687            keep <- character()
688            if(iLU <- names(yf) == "LU") {
689                ## TODO keep <- "LU"
690            }
691            y@factors <- yf[keep]
692        }
693        y
694    }
695    
696    setMethod("crossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
697              function(x, y = NULL) diagCspprod(x, y))
698    
699  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
700            function(x, y = NULL) crossprod(as(x, "sparseMatrix"), y))            function(x, y = NULL) diagCspprod(x, as(y, "CsparseMatrix")))
701    
702    ## Prefer calling diagCspprod to Cspdiagprod if going to transpose anyway
703    ##  x'y == (y'x)'
704    setMethod("crossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
705              function(x, y = NULL) t(diagCspprod(y, x)))
706    
707  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
708            function(x, y = NULL) crossprod(x, as(y, "sparseMatrix")))            function(x, y = NULL) t(diagCspprod(y, as(x, "Csparsematrix"))))
709    
710    setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
711              function(x, y = NULL) diagCspprod(x, t(y)))
712    
713  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
714            function(x, y = NULL) tcrossprod(as(x, "sparseMatrix"), y))            function(x, y = NULL) diagCspprod(x, t(as(y, "CsparseMatrix"))))
715    
716    setMethod("tcrossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
717              function(x, y = NULL) Cspdiagprod(x, y))
718    
719  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
720            function(x, y = NULL) tcrossprod(x, as(y, "sparseMatrix")))            function(x, y = NULL) Cspdiagprod(as(x, "CsparseMatrix"), y))
721    
722    setMethod("%*%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
723              function(x, y) diagCspprod(x, y))
724    
 ## FIXME?: In theory, this can be done *FASTER*, in some cases, via tapply1()  
725  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
726            function(x, y) as(x, "sparseMatrix") %*% y)            function(x, y) diagCspprod(as(x, "CsparseMatrix"), y))
727    
728  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
729            function(x, y) x %*% as(y, "sparseMatrix"))            function(x, y) Cspdiagprod(as(x, "CsparseMatrix"), y))
730  ## 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)  
731  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
732            function(x, y) x %*% as(y, "CsparseMatrix"))            function(x, y) Cspdiagprod(x, y))
733  ## NB: this is *not* needed for Tsparse & Rsparse  
734  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*
735  ##       do indeed work by going through sparse (and *not* ddense)!  ##       do indeed work by going through sparse (and *not* ddense)!
736    
   
   
737  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),
738            function(a, b, ...) {            function(a, b, ...) {
739                a@x <- 1/ a@x                a@x <- 1/ a@x
# Line 507  Line 742 
742            })            })
743    
744  solveDiag <- function(a, b, ...) {  solveDiag <- function(a, b, ...) {
745      if((n <- a@Dim[1]) != nrow(b))      if(a@Dim[1] != nrow(b))
746          stop("incompatible matrix dimensions")          stop("incompatible matrix dimensions")
747      ## trivially invert a 'in place' and multiply:      ## trivially invert a 'in place' and multiply:
748      a@x <- 1/ a@x      a@x <- 1/ a@x
# Line 526  Line 761 
761  ###---------------- <Ops> (<Arith>, <Logic>, <Compare> ) ----------------------  ###---------------- <Ops> (<Arith>, <Logic>, <Compare> ) ----------------------
762    
763  ## 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.)  
764  diagOdiag <- function(e1,e2) {  diagOdiag <- function(e1,e2) {
765      ## result should also be diagonal _ if possible _      ## result should also be diagonal _ if possible _
766      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 534  Line 768 
768      r00 <- callGeneric(if(is.numeric(e1@x)) 0 else FALSE,      r00 <- callGeneric(if(is.numeric(e1@x)) 0 else FALSE,
769                         if(is.numeric(e2@x)) 0 else FALSE)                         if(is.numeric(e2@x)) 0 else FALSE)
770      if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* diagonal      if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* diagonal
771          if(is.numeric(r)) {          if(is.numeric(r)) { # "double" *or* "integer"
772              if(is.numeric(e2@x)) {              if(is.numeric(e2@x)) {
773                  e2@x <- r; return(.diag.2N(e2)) }                  e2@x <- r; return(.diag.2N(e2)) }
774              if(!is.numeric(e1@x))              if(!is.numeric(e1@x))
775                  ## e.g. e1, e2 are logical;                  ## e.g. e1, e2 are logical;
776                  e1 <- as(e1, "dMatrix")                  e1 <- as(e1, "dMatrix")
777                if(!is.double(r)) r <- as.double(r)
778          }          }
779          else if(is.logical(r))          else if(is.logical(r))
780              e1 <- as(e1, "lMatrix")              e1 <- as(e1, "lMatrix")
781          else stop("intermediate 'r' is of type", typeof(r))          else stop(gettextf("intermediate 'r' is of type %s",
782                               typeof(r)), domain=NA)
783          e1@x <- r          e1@x <- r
784          .diag.2N(e1)          .diag.2N(e1)
785      }      }
786      else { ## result not diagonal, but at least symmetric:      else { ## result not diagonal, but at least symmetric:
787            ## e.g., m == m
788          isNum <- (is.numeric(r) || is.numeric(r00))          isNum <- (is.numeric(r) || is.numeric(r00))
789          isLog <- (is.logical(r) || is.logical(r00))          isLog <- (is.logical(r) || is.logical(r00))
790            Matrix.msg("exploding <diag> o <diag> into dense matrix", .M.level = 2)
         if(getOption("verbose"))  
             message("exploding  <diag>  o  <diag>  into dense matrix")  
791          d <- e1@Dim          d <- e1@Dim
792          n <- d[1]          n <- d[1]
793          stopifnot(length(r) == n)          stopifnot(length(r) == n)
794            if(isNum && !is.double(r)) r <- as.double(r)
795          xx <- as.vector(matrix(rbind(r, matrix(r00,n,n)), n,n))          xx <- as.vector(matrix(rbind(r, matrix(r00,n,n)), n,n))
796          newcl <-          newcl <-
797              paste(if(isNum) "d" else if(isLog) {              paste0(if(isNum) "d" else if(isLog) {
798                  if(!any(is.na(r)) && !any(is.na(r00))) "n" else "l"                  if(!any(is.na(r)) && !any(is.na(r00))) "n" else "l"
799              } else stop("not yet implemented .. please report")              } else stop("not yet implemented .. please report"), "syMatrix")
                   ,  
                   "syMatrix", sep='')  
800    
801          new(newcl, Dim = e1@Dim, Dimnames = e1@Dimnames, x = xx)          new(newcl, Dim = e1@Dim, Dimnames = e1@Dimnames, x = xx)
802      }      }
# Line 580  Line 814 
814              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)
815  }  }
816    
817  ## FIXME:    diagonal  o  triangular  |-->  triangular  ## diagonal  o  triangular  |-->  triangular
818  ## -----     diagonal  o  symmetric   |-->  symmetric  ## diagonal  o  symmetric   |-->  symmetric
819  ##    {also when other is sparse: do these "here" --  ##    {also when other is sparse: do these "here" --
820  ##     before conversion to sparse, since that loses "diagonality"}  ##     before conversion to sparse, since that loses "diagonality"}
821    diagOtri <- function(e1,e2) {
822        ## result must be triangular
823        r <- callGeneric(d1 <- .diag.x(e1), diag(e2)) # error if not "compatible"
824        ## Check what happens with non-diagonals, i.e. (0 o 0), (FALSE o 0), ...:
825        e1.0 <- if(.n1 <- is.numeric(d1   )) 0 else FALSE
826        r00 <- callGeneric(e1.0, if(.n2 <- is.numeric(e2[0])) 0 else FALSE)
827        if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* triangular
828            diag(e2) <- r
829            ## check what happens "in the triangle"
830            e2.2 <- if(.n2) 2 else TRUE
831            if(!callGeneric(e1.0, e2.2) == e2.2) { # values "in triangle" can change:
832                n <- dim(e2)[1L]
833                it <- indTri(n, upper = (e2@uplo == "U"))
834                e2[it] <- callGeneric(e1.0, e2[it])
835            }
836            e2
837        }
838        else { ## result not triangular ---> general
839            rr <- as(e2, "generalMatrix")
840            diag(rr) <- r
841            rr
842        }
843    }
844    
845    
846    setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "triangularMatrix"),
847              diagOtri)
848    ## For the reverse,  Ops == "Arith" | "Compare" | "Logic"
849    ##   'Arith'  :=  '"+"', '"-"', '"*"', '"^"', '"%%"', '"%/%"', '"/"'
850    setMethod("Arith", signature(e1 = "triangularMatrix", e2 = "diagonalMatrix"),
851              function(e1,e2)
852          { ## this must only trigger for *dense* e1
853              switch(.Generic,
854                     "+" = .Call(dtrMatrix_addDiag, as(e1,"dtrMatrix"),   .diag.x(e2)),
855                     "-" = .Call(dtrMatrix_addDiag, as(e1,"dtrMatrix"), - .diag.x(e2)),
856                     "*" = {
857                         n <- e2@Dim[1L]
858                         d2 <- if(e2@diag == "U") { # unit-diagonal
859                             d <- rep.int(as1(e2@x), n)
860                             e2@x <- d
861                             e2@diag <- "N"
862                             d
863                         } else e2@x
864                         e2@x <- diag(e1) * d2
865                         e2
866                     },
867                     "^" = { ## will be dense ( as  <ANY> ^ 0 == 1 ):
868                         e1 ^ as(e2, "denseMatrix")
869                     },
870                     ## otherwise:
871                     callGeneric(e1, diag2Tsmart(e2,e1)))
872    })
873    
874    ## Compare --> 'swap' (e.g.   e1 < e2   <==>  e2 > e1 ):
875    setMethod("Compare", signature(e1 = "triangularMatrix", e2 = "diagonalMatrix"),
876              .Cmp.swap)
877    ## '&' and "|'  are commutative:
878    setMethod("Logic", signature(e1 = "triangularMatrix", e2 = "diagonalMatrix"),
879              function(e1,e2) callGeneric(e2, e1))
880    
881  ## For almost everything else, diag* shall be treated "as sparse" :  ## For almost everything else, diag* shall be treated "as sparse" :
882  ## These are cheap implementations via coercion  ## These are cheap implementations via coercion
# Line 607  Line 900 
900  ##       Compare --> logical  ##       Compare --> logical
901  ##       Logic   --> logical: "lMatrix"  ##       Logic   --> logical: "lMatrix"
902    
903  ##  other = "numeric" : stay diagonal if possible  ## Other = "numeric" : stay diagonal if possible
904  ## ddi*: Arith: result numeric, potentially ddiMatrix  ## ddi*: Arith: result numeric, potentially ddiMatrix
905  setMethod("Arith", signature(e1 = "ddiMatrix", e2 = "numeric"),  for(arg2 in c("numeric","logical"))
906    setMethod("Arith", signature(e1 = "ddiMatrix", e2 = arg2),
907            function(e1,e2) {            function(e1,e2) {
908                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]
909                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
910                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
911                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
912                    if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)                    if(e1@diag == "U") {
913                    if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {                        if(any((r <- callGeneric(1, e2)) != 1)) {
914                        e1@diag <- "N"                        e1@diag <- "N"
915                        if(L1) r <- rep.int(r, n)                            e1@x[seq_len(n)] <- r # possibly recycling r
916                    } else                        } ## else: result = e1  (is "U" diag)
917                      } else {
918                        r <- callGeneric(e1@x, e2)                        r <- callGeneric(e1@x, e2)
919                    e1@x <- if(L1) r else r[1L + n*(0:(n-1L))]                        ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
920                    return(e1)                        e1@x[] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
921                }                }
922                      e1
923                  } else
924                callGeneric(diag2tT.u(e1,e2, "d"), e2)                callGeneric(diag2tT.u(e1,e2, "d"), e2)
925            })            })
926    
927  setMethod("Arith", signature(e1 = "numeric", e2 = "ddiMatrix"),  for(arg1 in c("numeric","logical"))
928    setMethod("Arith", signature(e1 = arg1, e2 = "ddiMatrix"),
929            function(e1,e2) {            function(e1,e2) {
930                n <- e2@Dim[1]; nsq <- n*n                n <- e2@Dim[1]
931                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
932                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
933                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
934                    if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)                    if(e2@diag == "U") {
935                    if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {                        if(any((r <- callGeneric(e1, 1)) != 1)) {
936                        e2@diag <- "N"                        e2@diag <- "N"
937                        if(L1) r <- rep.int(r, n)                            e2@x[seq_len(n)] <- r # possibly recycling r
938                    } else                        } ## else: result = e2  (is "U" diag)
939                      } else {
940                        r <- callGeneric(e1, e2@x)                        r <- callGeneric(e1, e2@x)
941                    e2@x <- if(L1) r else r[1L + n*(0:(n-1L))]                        ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
942                    return(e2)                        e2@x[] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
943                }                }
944                      e2
945                  } else
946                callGeneric(e1, diag2tT.u(e2,e1, "d"))                callGeneric(e1, diag2tT.u(e2,e1, "d"))
947            })            })
948    
949  ## ldi* Arith --> result numeric, potentially ddiMatrix  ## ldi* Arith --> result numeric, potentially ddiMatrix
950  setMethod("Arith", signature(e1 = "ldiMatrix", e2 = "numeric"),  for(arg2 in c("numeric","logical"))
951    setMethod("Arith", signature(e1 = "ldiMatrix", e2 = arg2),
952            function(e1,e2) {            function(e1,e2) {
953                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]
954                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
955                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
956                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
957                    if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)                    E <- copyClass(e1, "ddiMatrix", c("diag", "Dim", "Dimnames"))#FIXME: if ok, check=FALSE
958                    if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {                    ## E <- copyClass(e1, "ddiMatrix", check=FALSE)
959                        e1@diag <- "N"                    ## storage.mode(E@x) <- "double"
960                        if(L1) r <- rep.int(r, n)                    if(e1@diag == "U") {
961                    } else                        if(any((r <- callGeneric(1, e2)) != 1)) {
962                              E@diag <- "N"
963                              E@x[seq_len(n)] <- r # possibly recycling r
964                          } ## else: result = E  (is "U" diag)
965                      } else {
966                        r <- callGeneric(e1@x, e2)                        r <- callGeneric(e1@x, e2)
967                    e1 <- copyClass(e1, "ddiMatrix", c("diag", "Dim", "Dimnames"))                        ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
968                    e1@x <- if(L1) r else r[1L + n*(0:(n-1L))]                        E@x[seq_len(n)] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
                   return(e1)  
969                }                }
970                callGeneric(diag2tT.u(e1,e2, "d"), e2)                    E
971                  } else
972                      callGeneric(diag2tT.u(e1,e2, "l"), e2)
973            })            })
974    
975  setMethod("Arith", signature(e1 = "numeric", e2 = "ldiMatrix"),  for(arg1 in c("numeric","logical"))
976    setMethod("Arith", signature(e1 = arg1, e2 = "ldiMatrix"),
977            function(e1,e2) {            function(e1,e2) {
978                n <- e2@Dim[1]; nsq <- n*n                n <- e2@Dim[1]
979                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
980                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
981                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
982                    if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)                    E <- copyClass(e2, "ddiMatrix", c("diag", "Dim", "Dimnames"))#FIXME: if ok, check=FALSE
983                    if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {                    ## E <- copyClass(e2, "ddiMatrix", check=FALSE)
984                        e2@diag <- "N"                    ## storage.mode(E@x) <- "double"
985                        if(L1) r <- rep.int(r, n)                    if(e2@diag == "U") {
986                    } else                        if(any((r <- callGeneric(e1, 1)) != 1)) {
987                              E@diag <- "N"
988                              E@x[seq_len(n)] <- r # possibly recycling r
989                          } ## else: result = E  (is "U" diag)
990                      } else {
991                        r <- callGeneric(e1, e2@x)                        r <- callGeneric(e1, e2@x)
992                    e2 <- copyClass(e2, "ddiMatrix", c("diag", "Dim", "Dimnames"))                        ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
993                    e2@x <- if(L1) r else r[1L + n*(0:(n-1L))]                        E@x[seq_len(n)] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
                   return(e2)  
994                }                }
995                callGeneric(e1, diag2tT.u(e2,e1, "d"))                    E
996                  } else
997                      callGeneric(e1, diag2tT.u(e2,e1, "l"))
998            })            })
999    
1000  ## ddi*: for "Ops" without Arith --> result logical, potentially ldi  ## ddi*: for "Ops" without "Arith": <Compare> or <Logic> --> result logical, potentially ldi
1001  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "numeric"),  ##
1002    ## Note that  ("numeric", "ddiMatrix")  is simply swapped, e.g.,
1003    if(FALSE) {
1004        selectMethod("<", c("numeric","lMatrix"))# Compare
1005        selectMethod("&", c("numeric","lMatrix"))# Logic
1006    } ## so we don't need to define a method here :
1007    
1008    for(arg2 in c("numeric","logical"))
1009    setMethod("Ops", signature(e1 = "ddiMatrix", e2 = arg2),
1010            function(e1,e2) {            function(e1,e2) {
1011                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]
1012                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
1013                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
1014                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
1015                    if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)                    E <- copyClass(e1, "ldiMatrix", c("diag", "Dim", "Dimnames"))#FIXME: if ok, check=FALSE
1016                    if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {                    ## E <- copyClass(e1, "ldiMatrix", check=FALSE)
1017                        e1@diag <- "N"                    ## storage.mode(E@x) <- "logical"
1018                        if(L1) r <- rep.int(r, n)                    if(e1@diag == "U") {
1019                    } else                        if(any((r <- callGeneric(1, e2)) != 1)) {
1020                              E@diag <- "N"
1021                              E@x[seq_len(n)] <- r # possibly recycling r
1022                          } ## else: result = E  (is "U" diag)
1023                      } else {
1024                        r <- callGeneric(e1@x, e2)                        r <- callGeneric(e1@x, e2)
1025                    e1 <- copyClass(e1, "ldiMatrix", c("diag", "Dim", "Dimnames"))                        ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
1026                    e1@x <- if(L1) r else r[1L + n*(0:(n-1L))]                        E@x[seq_len(n)] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
                   return(e1)  
1027                }                }
1028                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)  
1029                    } else                    } else
1030                        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*(0:(n-1L))]  
                   return(e2)  
               }  
               callGeneric(e1, diag2tT.u(e2,e1, "l"))  
1031            })            })
1032    
1033  ## ldi*: for "Ops" without Arith --> result logical, potentially ldi  ## ldi*: for "Ops" without "Arith": <Compare> or <Logic> --> result logical, potentially ldi
1034  setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "numeric"),  for(arg2 in c("numeric","logical"))
1035    setMethod("Ops", signature(e1 = "ldiMatrix", e2 = arg2),
1036            function(e1,e2) {            function(e1,e2) {
1037                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]
1038                f0 <- callGeneric(FALSE, e2)                f0 <- callGeneric(FALSE, e2)
1039                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
1040                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
1041                    if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)  
1042                    if(e1@diag == "U" && any((r <- callGeneric(TRUE, e2)) != 1)) {                    if(e1@diag == "U") {
1043                          if(any((r <- callGeneric(TRUE, e2)) != 1)) {
1044                        e1@diag <- "N"                        e1@diag <- "N"
1045                        if(L1) r <- rep.int(r, n)                            e1@x[seq_len(n)] <- r # possibly recycling r
1046                    } else                        } ## else: result = e1  (is "U" diag)
1047                      } else {
1048                        r <- callGeneric(e1@x, e2)                        r <- callGeneric(e1@x, e2)
1049                    e1@x <- if(L1) r else r[1L + n*(0:(n-1L))]                        ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
1050                    return(e1)                        e1@x[] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
1051                }                }
1052                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)  
1053                    } else                    } else
1054                        r <- callGeneric(e1, e2@x)                    callGeneric(diag2tT.u(e1,e2, "l"), e2)
                   e2@x <- if(L1) r else r[1L + n*(0:(n-1L))]  
                   return(e2)  
               }  
               callGeneric(e1, diag2tT.u(e2,e1, "l"))  
1055            })            })
1056    
1057    
   
1058  ## Not {"sparseMatrix", "numeric} :  {"denseMatrix", "matrix", ... }  ## Not {"sparseMatrix", "numeric} :  {"denseMatrix", "matrix", ... }
1059    for(other in c("ANY", "Matrix", "dMatrix")) {
1060  ## ddi*:  ## ddi*:
1061  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "ANY"),      setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other),
1062            function(e1,e2) callGeneric(diag2tT.u(e1,e2, "d"), e2))                function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2))
1063  setMethod("Ops", signature(e1 = "ANY", e2 = "ddiMatrix"),      setMethod("Ops", signature(e1 = other, e2 = "ddiMatrix"),
1064            function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "d")))                function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "d")))
1065  ## ldi*:  ## ldi*:
1066  setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "ANY"),      setMethod("Ops", signature(e1 = "ldiMatrix", e2 = other),
1067            function(e1,e2) callGeneric(diag2tT.u(e1,e2, "l"), e2))                function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2))
1068  setMethod("Ops", signature(e1 = "ANY", e2 = "ldiMatrix"),      setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),
1069            function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "l")))                function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "l")))
1070    }
1071    
1072  ## This should *not* dispatch to <dense> methods (in ./Ops.R ), as  ## Direct subclasses of "denseMatrix": currently ddenseMatrix, ldense... :
1073  ##  FALSE & <anything> |-> FALSE : hence result should be diagonal:  dense.subCl <- local({ dM.scl <- getClass("denseMatrix")@subclasses
1074  for(cl in diCls) {                         names(dM.scl)[vapply(dM.scl, slot, 0, "distance") == 1] })
1075      setMethod("&", signature(e1 = cl, e2 = "ANY"),  for(DI in diCls) {
1076                function(e1,e2) e1 & as(e2,"Matrix"))      dMeth <- if(extends(DI, "dMatrix"))
1077      setMethod("&", signature(e1 = "ANY", e2 = cl),          function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2)
1078                function(e1,e2) as(e1,"Matrix") & e2)      else # "lMatrix", the only other kind for now
1079      for(c2 in c("denseMatrix", "Matrix")) {          function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2)
1080          setMethod("&", signature(e1 = cl, e2 = c2),      for(c2 in c(dense.subCl, "Matrix")) {
1081                    function(e1,e2) e1 & Diagonal(x = diag(e2)))          for(Fun in c("*", "&")) {
1082          setMethod("&", signature(e1 = c2, e2 = cl),              setMethod(Fun, signature(e1 = DI, e2 = c2),
1083                    function(e1,e2) Diagonal(x = diag(e1)) & e2)                        function(e1,e2) callGeneric(e1, Diagonal(x = diag(e2))))
1084                setMethod(Fun, signature(e1 = c2, e2 = DI),
1085                          function(e1,e2) callGeneric(Diagonal(x = diag(e1)), e2))
1086            }
1087            setMethod("^", signature(e1 = c2, e2 = DI),
1088                      function(e1,e2) callGeneric(Diagonal(x = diag(e1)), e2))
1089            for(Fun in c("%%", "%/%", "/")) ## 0 <op> 0 |--> NaN  for these.
1090                setMethod(Fun, signature(e1 = DI, e2 = c2), dMeth)
1091      }      }
1092  }  }
1093    
1094    
1095  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"
1096  ### ----------  any, all: separately here  ### ----------   the last 4: separately here
1097  for(cl in diCls) {  for(cl in diCls) {
1098  setMethod("any", cl,  setMethod("any", cl,
1099            function (x, ..., na.rm) {            function (x, ..., na.rm) {
1100                if(any(x@Dim == 0)) FALSE                if(any(x@Dim == 0)) FALSE
1101                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)
1102            })            })
1103  setMethod("all",  cl, function (x, ..., na.rm) any(x@Dim == 0))  setMethod("all",  cl, function (x, ..., na.rm) {
1104  setMethod("prod", cl, function (x, ..., na.rm) as.numeric(any(x@Dim == 0)))      n <- x@Dim[1]
1105        if(n >= 2) FALSE
1106        else if(n == 0 || x@diag == "U") TRUE
1107        else all(x@x, ..., na.rm = na.rm)
1108    })
1109    setMethod("prod", cl, function (x, ..., na.rm) {
1110        n <- x@Dim[1]
1111        if(n >= 2) 0
1112        else if(n == 0 || x@diag == "U") 1
1113        else ## n == 1, diag = "N" :
1114            prod(x@x, ..., na.rm = na.rm)
1115    })
1116    
1117  setMethod("sum", cl,  setMethod("sum", cl,
1118            function(x, ..., na.rm) {            function(x, ..., na.rm) {
# Line 838  Line 1151 
1151      invisible(x)      invisible(x)
1152  }  }
1153    
1154    ## somewhat consistent with "print" for sparseMatrix :
1155    setMethod("print", signature(x = "diagonalMatrix"), prDiag)
1156    
1157  setMethod("show", signature(object = "diagonalMatrix"),  setMethod("show", signature(object = "diagonalMatrix"),
1158            function(object) {            function(object) {
1159                d <- dim(object)                d <- dim(object)
1160                cl <- class(object)                cl <- class(object)
1161                cat(sprintf('%d x %d diagonal matrix of class "%s"\n',                cat(sprintf('%d x %d diagonal matrix of class "%s"',
1162                            d[1], d[2], cl))                            d[1], d[2], cl))
1163                  if(d[1] < 50) {
1164                      cat("\n")
1165                prDiag(object)                prDiag(object)
1166                  } else {
1167                      cat(", with diagonal entries\n")
1168                      show(diag(object))
1169                      invisible(object)
1170                  }
1171            })            })
1172    
1173    rm(arg1, arg2, other, DI, cl, c1, c2,
1174       dense.subCl, diCls)# not used elsewhere
1175    
1176    setMethod("summary", signature(object = "diagonalMatrix"),
1177              function(object, ...) {
1178                  d <- dim(object)
1179                  r <- summary(object@x, ...)
1180                  attr(r, "header") <-
1181                      sprintf('%d x %d diagonal Matrix of class "%s"',
1182                              d[1], d[2], class(object))
1183                  ## use ole' S3 technology for such a simple case
1184                  class(r) <- c("diagSummary", class(r))
1185                  r
1186              })
1187    
1188    print.diagSummary <- function (x, ...) {
1189        cat(attr(x, "header"),"\n")
1190        class(x) <- class(x)[-1]
1191        print(x, ...)
1192        invisible(x)
1193    }

Legend:
Removed from v.2185  
changed lines
  Added in v.2904

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