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 1805, Tue Mar 27 16:46:03 2007 UTC pkg/Matrix/R/diagMatrix.R revision 2984, Sat Apr 12 21:37:37 2014 UTC
# Line 5  Line 5 
5  ##          but *not* diag() extractor!  ##          but *not* diag() extractor!
6  Diagonal <- function(n, x = NULL)  Diagonal <- function(n, x = NULL)
7  {  {
8      ## Allow  Diagonal(4)  and  Diagonal(x=1:5)      ## Allow  Diagonal(4), Diagonal(x=1:5), and  Diagonal(4, TRUE)
9      if(missing(n))      n <- if(missing(n)) length(x) else {
         n <- length(x)  
     else {  
10          stopifnot(length(n) == 1, n == as.integer(n), n >= 0)          stopifnot(length(n) == 1, n == as.integer(n), n >= 0)
11          n <- as.integer(n)          as.integer(n)
12      }      }
13    
14      if(missing(x)) ## unit diagonal matrix      if(missing(x)) ## unit diagonal matrix
15          new("ddiMatrix", Dim = c(n,n), diag = "U")          new("ddiMatrix", Dim = c(n,n), diag = "U")
16      else {      else {
17          stopifnot(length(x) == n)          lx <- length(x)
18            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 26  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          new(cl, Dim = c(n,n), diag = "N", x = x)          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",
33                    x = if(lx.1) rep.int(x,n) else x)
34      }      }
35  }  }
36    
37  ### This is modified from a post of Bert Gunter to R-help on  1 Sep 2005.  .sparseDiagonal <- function(n, x = 1, uplo = "U",
38  ### Bert's code built on a post by Andy Liaw who most probably was influenced                              shape = if(missing(cols)) "t" else "g",
39  ### by earlier posts, notably one by Scott Chasalow on S-news, 16 Jan 2002                              unitri, kind,
40  ### who posted his bdiag() function written in December 1995.                              cols = if(n) 0:(n - 1L) else integer(0))
41    {
42        stopifnot(n == (n. <- as.integer(n)), (n <- n.) >= 0)
43        if(!(mcols <- missing(cols)))
44            stopifnot(0 <= (cols <- as.integer(cols)), cols < n)
45        m <- length(cols)
46        if(missing(kind))
47            kind <-
48                if(is.double(x)) "d"
49                else if(is.logical(x)) "l"
50                else { ## for now
51                    storage.mode(x) <- "double"
52                    "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        }
79    }
80    
81  bdiag <- function(...) {  ## Pkg 'spdep' had (relatively slow) versions of this as_dsCMatrix_I()
82      if(nargs() == 0) return(new("dgCMatrix"))  .symDiagonal <- function(n, x = rep.int(1,n), uplo = "U")
83      ## else :      .sparseDiagonal(n, x, uplo, shape = "s")
84      mlist <- if (nargs() == 1) as.list(...) else list(...)  
85      dims <- sapply(mlist, dim)  # instead of   diagU2N(as(Diagonal(n), "CsparseMatrix")), diag = "N" in any case:
86    .trDiagonal <- function(n, x = 1, uplo = "U", unitri=TRUE)
87        .sparseDiagonal(n, x, uplo, shape = "t", unitri=unitri)
88    
89    
90    ## This is modified from a post of Bert Gunter to R-help on  1 Sep 2005.
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(0:0, 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, 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  diag2T <- function(from) {  bdiag <- function(...) {
190      i <- if(from@diag == "U") integer(0) else seq_len(from@Dim[1]) - 1:1      if((nA <- nargs()) == 0) return(new("dgCMatrix"))
191      new(paste(.M.kind(from), "tTMatrix", sep=''),      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)) {
203        ## to triangular Tsparse
204        i <- if(from@diag == "U") integer(0) else seq_len(from@Dim[1]) - 1L
205        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,
208          x = from@x, # <- ok for diag = "U" and "N" (!)          x = from@x, # <- ok for diag = "U" and "N" (!)
209          i = i, j = i)          i = i, j = i)
210  }  }
211    
212  setAs("diagonalMatrix", "triangularMatrix", diag2T)  .diag2sT <- function(from, uplo = "U", kind = .M.kind(from)) {
213  setAs("diagonalMatrix", "sparseMatrix", diag2T)      ## to symmetric Tsparse
214        n <- from@Dim[1]
215        i <- seq_len(n) - 1L
216        new(paste0(kind, "sTMatrix"),
217            Dim = from@Dim, Dimnames = from@Dimnames,
218            i = i, j = i, uplo = uplo,
219            x = if(from@diag == "N") from@x else ## "U"-diag
220            rep.int(switch(kind,
221                           "d" = 1.,
222                           "l" =,
223                           "n" = TRUE,
224                           ## otherwise
225                           stop(gettextf("%s kind not yet implemented",
226                                         sQuote(kind)), domain=NA)),
227                    n))
228    }
229    
230    ## diagonal -> triangular,  upper / lower depending on "partner":
231    diag2tT.u <- function(d, x, kind = .M.kind(d))
232        .diag2tT(d, uplo = if(is(x,"triangularMatrix")) x@uplo else "U", kind)
233    
234    ## diagonal -> sparse {triangular OR symmetric} (upper / lower) depending on "partner":
235    diag2Tsmart <- function(d, x, kind = .M.kind(d)) {
236        clx <- getClassDef(class(x))
237        if(extends(clx, "symmetricMatrix"))
238            .diag2sT(d, uplo = x@uplo, kind)
239        else
240            .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,
250    ## and because we can save a .M.kind() call, we use this explicit
251    ## "hack"  instead of signature  x = "diagonalMatrix" :
252    ##
253    ## ddi*:
254    di2tT <- function(from) .diag2tT(from, "U", "d")
255    setAs("ddiMatrix", "triangularMatrix", di2tT)
256    ##_no_longer_ setAs("ddiMatrix", "sparseMatrix", di2tT)
257  ## needed too (otherwise <dense> -> Tsparse is taken):  ## needed too (otherwise <dense> -> Tsparse is taken):
258  setAs("diagonalMatrix", "TsparseMatrix", diag2T)  setAs("ddiMatrix", "TsparseMatrix", di2tT)
259  ## is better than this:  setAs("ddiMatrix", "dsparseMatrix", di2tT)
260  ## setAs("diagonalMatrix", "sparseMatrix",  setAs("ddiMatrix", "CsparseMatrix",
261  ##       function(from)        function(from) as(.diag2tT(from, "U", "d"), "CsparseMatrix"))
262  ##        as(from, if(is(from, "dMatrix")) "dgCMatrix" else "lgCMatrix"))  setAs("ddiMatrix", "symmetricMatrix",
263  setAs("diagonalMatrix", "CsparseMatrix",        function(from) .diag2sT(from, "U", "d"))
264        function(from) as(diag2T(from), "CsparseMatrix"))  ##
265    ## ldi*:
266    di2tT <- function(from) .diag2tT(from, "U", "l")
267    setAs("ldiMatrix", "triangularMatrix", di2tT)
268    ##_no_longer_ setAs("ldiMatrix", "sparseMatrix", di2tT)
269    ## needed too (otherwise <dense> -> Tsparse is taken):
270    setAs("ldiMatrix", "TsparseMatrix", di2tT)
271    setAs("ldiMatrix", "lsparseMatrix", di2tT)
272    setAs("ldiMatrix", "CsparseMatrix",
273          function(from) as(.diag2tT(from, "U", "l"), "CsparseMatrix"))
274    setAs("ldiMatrix", "symmetricMatrix",
275          function(from) .diag2sT(from, "U", "l"))
276    rm(di2tT)
277    
278  setAs("diagonalMatrix", "matrix",  setAs("diagonalMatrix", "nMatrix",
279        function(from) {        function(from) {
280            n <- from@Dim[1]            n <- from@Dim[1]
281            diag(x = if(from@diag == "U") { if(is.numeric(from@x)) 1. else TRUE            i <- if(from@diag == "U") integer(0) else which(isN0(from@x)) - 1L
282                                       } else from@x,            new("ntTMatrix", i = i, j = i, diag = from@diag,
283                 nrow = n, ncol = n)                Dim = from@Dim, Dimnames = from@Dimnames)
284        })        })
285    
286  setAs("diagonalMatrix", "generalMatrix", # prefer sparse:  setAs("diagonalMatrix", "nsparseMatrix", function(from) as(from, "nMatrix"))
       function(from) as(from, paste(.M.kind(from), "gCMatrix", sep='')))  
   
 .diag.x <- function(m) {  
     if(m@diag == "U")  
         rep.int(if(is.numeric(m@x)) 1. else TRUE,  
                 m@Dim[1])  
     else m@x  
 }  
287    
288  .diag.2N <- function(m) {  ## Cheap fast substitute for diag() which *does* preserve the mode of x :
289      if(m@diag == "U") m@diag <- "N"  mkDiag <- function(x, n) {
290      m      y <- matrix(as0(mod=mode(x)), n,n)
291        if (n > 0) y[1L + 0:(n - 1L) * (n + 1)] <- x
292        y
293  }  }
294    
295  ## given the above, the following  4  coercions should be all unneeded;  setAs("diagonalMatrix", "matrix",
 ## we prefer triangular to general:  
 setAs("ddiMatrix", "dgTMatrix",  
296        function(from) {        function(from) {
297            .Deprecated("as(, \"sparseMatrix\")")            ## want "ldiMatrix" -> <logical> "matrix" :
298            n <- from@Dim[1]            mkDiag(if(from@diag == "U") as1(from@x) else from@x,
299            i <- seq_len(n) - 1:1                   n = from@Dim[1])
300            new("dgTMatrix", i = i, j = i, x = .diag.x(from),        })
               Dim = c(n,n), Dimnames = from@Dimnames) })  
301    
302  setAs("ddiMatrix", "dgCMatrix",  setMethod("as.vector", signature(x = "diagonalMatrix", mode="missing"),
303        function(from) as(as(from, "sparseMatrix"), "dgCMatrix"))            function(x, mode) {
304                  n <- x@Dim[1]
305                  mod.x <- mode(x@x)
306                  r <- vector(mod.x, length = n^2)
307                  if(n)
308                      r[1 + 0:(n - 1L) * (n + 1)] <-
309                          if(x@diag == "U") as1(mod=mod.x) else x@x
310                  r
311              })
312    
313  setAs("ldiMatrix", "lgTMatrix",  setAs("diagonalMatrix", "generalMatrix", # prefer sparse:
314        function(from) {        function(from) as(as(from, "CsparseMatrix"), "generalMatrix"))
           .Deprecated("as(, \"sparseMatrix\")")  
           n <- from@Dim[1]  
           if(from@diag == "U") { # unit-diagonal  
               x <- rep.int(TRUE, n)  
               i <- seq_len(n) - 1:1  
           } else { # "normal"  
               nz <- nz.NA(from@x, na. = TRUE)  
               x <- from@x[nz]  
               i <- which(nz) - 1:1  
           }  
           new("lgTMatrix", i = i, j = i, x = x,  
               Dim = c(n,n), Dimnames = from@Dimnames) })  
315    
316  setAs("ldiMatrix", "lgCMatrix",  setAs("diagonalMatrix", "denseMatrix",
317        function(from) as(as(from, "lgTMatrix"), "lgCMatrix"))        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
321    
322    .diag.2N <- function(m) {
323        if(m@diag == "U") m@diag <- "N"
324        m
325    }
326    
 if(FALSE) # now have faster  "ddense" -> "dge"  
327  setAs("ddiMatrix", "dgeMatrix",  setAs("ddiMatrix", "dgeMatrix",
328        function(from) as(as(from, "matrix"), "dgeMatrix"))        function(from) .Call(dup_mMatrix_as_dgeMatrix, from))
329    setAs("ddiMatrix", "ddenseMatrix",
330          function(from) as(as(from, "triangularMatrix"),"denseMatrix"))
331    setAs("ldiMatrix", "ldenseMatrix",
332          function(from) as(as(from, "triangularMatrix"),"denseMatrix"))
333    
334    
335  setAs("matrix", "diagonalMatrix",  setAs("matrix", "diagonalMatrix",
336        function(from) {        function(from) {
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 165  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        })        })
373    
374    
375  setMethod("diag", signature(x = "diagonalMatrix"),  setMethod("diag", signature(x = "diagonalMatrix"),
376            function(x = 1, nrow, ncol = n) .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(isDiagonal(x)) as(x, "diagonalMatrix") else x      if(isS4(x) && isDiagonal(x)) as(x, "diagonalMatrix") else x
387  }  }
388    
389  setMethod("[", signature(x = "diagonalMatrix", i = "index",  setMethod("[", signature(x = "diagonalMatrix", i = "index",
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, 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, j, drop) subDiag(x, j=j, drop=drop))            function(x, i, j, ..., drop) subDiag(x, j=j, drop=drop))
403    
404  ## When you assign to a diagonalMatrix, the result should be  ## When you assign to a diagonalMatrix, the result should be
405  ## diagonal or sparse ---  ## diagonal or sparse ---
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  replDiag <- function(x, i, j, value) {  ## Only(?) current bug:  x[i] <- value  is wrong when  i is *vector*
408      x <- as(x, "sparseMatrix")  replDiag <- function(x, i, j, ..., value) {
409        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))      else if(missing(j)) { ##  x[i , ] <- v  *OR*   x[i] <- v
413            na <- nargs()
414    ##         message("diagnosing replDiag() -- nargs()= ", na)
415            if(na == 4)
416          x[i, ] <- value          x[i, ] <- value
417      else          else if(na == 3)
418                x[i] <- value
419            else stop(gettextf("Internal bug: nargs()=%d; please report",
420                               na), domain=NA)
421        } 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
424  }  }
425    
426  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index",
427                                  j = "index", value = "replValue"), replDiag)                                  j = "index", value = "replValue"), replDiag)
428    
429  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index",
430                                  j = "missing", value = "replValue"),                                  j = "missing", value = "replValue"),
431                   function(x, i, value) replDiag(x, i=i, value=value))                   function(x,i,j, ..., value) {
432                         ## message("before replDiag() -- nargs()= ", nargs())
433                         if(nargs() == 3)
434                             replDiag(x, i=i, value=value)
435                         else ## nargs() == 4 :
436                             replDiag(x, i=i, , value=value)
437                     })
438    
439  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing",
440                                  j = "index", value = "replValue"),                                  j = "index", value = "replValue"),
441                   function(x, j, value) replDiag(x, j=j, value=value))                   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"),
463                     function(x,i,j, ..., value) {
464                         if(ncol(i) == 2) {
465                             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
474                                 x
475                             } else { ## no longer diagonal, but remain sparse:
476                                 x <- as(x, "TsparseMatrix")
477                                 x[i] <- value
478                                 x
479                             }
480                         }
481                         else { # behave as "base R": use as if vector
482                             x <- as(x, "matrix")
483                             x[i] <- value
484                             Matrix(x)
485                         }
486                     })
487    
488    
489    ## value = "sparseMatrix":
490    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)
521            function(object) TRUE)  
522  setMethod("isSymmetric", signature(object = "diagonalMatrix"),  setMethod("symmpart", signature(x = "diagonalMatrix"), function(x) x)
523            function(object) TRUE)  setMethod("skewpart", signature(x = "diagonalMatrix"), setZero)
524    
525  setMethod("chol", signature(x = "ddiMatrix"),# pivot = "ANY"  setMethod("chol", signature(x = "ddiMatrix"),
526            function(x, pivot) {            function(x, pivot, ...) {
527                if(any(x@x < 0)) stop("chol() is undefined for diagonal matrix with negative entries")                if(x@diag == "U") return(x)
528                  ## else
529                  if(any(x@x < 0))
530                      stop("chol() is undefined for diagonal matrix with negative entries")
531                x@x <- sqrt(x@x)                x@x <- sqrt(x@x)
532                x                x
533            })            })
534  ## chol(L) is L for logical diagonal:  ## chol(L) is L for logical diagonal:
535  setMethod("chol", signature(x = "ldiMatrix"), function(x, pivot) x)  setMethod("chol", signature(x = "ldiMatrix"), function(x, pivot, ...) x)
536    
537  setMethod("!", "ldiMatrix", function(e1) {  setMethod("determinant", signature(x = "diagonalMatrix", logarithm = "logical"),
538      if(e1@diag == "N")            function(x, logarithm, ...) mkDet(.diag.x(x), logarithm))
539          e1@x <- !e1@x  
540      else { ## "U"  setMethod("norm", signature(x = "diagonalMatrix", type = "character"),
541          e1@diag <- "N"            function(x, type, ...) {
542          e1@x <- rep.int(FALSE, e1@Dim[1])                if((n <- x@Dim[1]) == 0) return(0) # as for "sparseMatrix"
543                  type <- toupper(substr(type[1], 1, 1))
544                  isU <- (x@diag == "U") # unit-diagonal
545                  if(type == "F") sqrt(if(isU) n else sum(x@x^2))
546                  else { ## norm == "I","1","O","M" :
547                      if(isU) 1 else max(abs(x@x))
548      }      }
     e1  
549  })  })
550    
551    
552    
553  ## Basic Matrix Multiplication {many more to add}  ## Basic Matrix Multiplication {many more to add}
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      if(any(dim(x) != dim(y))) stop("non-matching dimensions")      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 285  Line 582 
582    
583    
584  diagmatprod <- function(x, y) {  diagmatprod <- function(x, y) {
585        ## x is diagonalMatrix
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"),
592            diagmatprod)            diagmatprod)
593    ## sneaky .. :
594  formals(diagmatprod) <- alist(x=, y=NULL)  formals(diagmatprod) <- alist(x=, y=NULL)
595  setMethod("crossprod", signature(x = "diagonalMatrix", y = "matrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "matrix"),
596            diagmatprod)            diagmatprod)
 setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "matrix"),  
           diagmatprod)  
597    
598  diagdgeprod <- function(x, y) {  diagGeprod <- function(x, y) {
599      dx <- dim(x)      dx <- dim(x)
600      dy <- dim(y)      dy <- dim(y)
601      if(dx[2] != dy[1]) stop("non-matching dimensions")      if(dx[2] != dy[1]) stop("non-matching dimensions")
# Line 308  Line 603 
603          y@x <- x@x * y@x          y@x <- x@x * y@x
604      y      y
605  }  }
606  setMethod("%*%", signature(x = "diagonalMatrix", y = "dgeMatrix"),  setMethod("%*%", signature(x= "diagonalMatrix", y= "dgeMatrix"), diagGeprod)
607            diagdgeprod, valueClass = "dgeMatrix")  setMethod("%*%", signature(x= "diagonalMatrix", y= "lgeMatrix"), diagGeprod)
608  formals(diagdgeprod) <- alist(x=, y=NULL)  formals(diagGeprod) <- alist(x=, y=NULL)
609  setMethod("crossprod", signature(x = "diagonalMatrix", y = "dgeMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "dgeMatrix"),
610            diagdgeprod, valueClass = "dgeMatrix")            diagGeprod, valueClass = "dgeMatrix")
611    setMethod("crossprod", signature(x = "diagonalMatrix", y = "lgeMatrix"),
612              diagGeprod)
613    
614  setMethod("%*%", signature(x = "matrix", y = "diagonalMatrix"),  matdiagprod <- function(x, y) {
           function(x, y) {  
615                dx <- dim(x)                dx <- dim(x)
616                dy <- dim(y)                dy <- dim(y)
617                if(dx[2] != dy[1]) stop("non-matching dimensions")                if(dx[2] != dy[1]) stop("non-matching dimensions")
618                as(if(y@diag == "U") x else x * rep(y@x, each = dx[1]), "Matrix")      Matrix(if(y@diag == "U") x else x * rep(y@x, each = dx[1]))
619            })  }
620    setMethod("%*%", signature(x = "matrix", y = "diagonalMatrix"),
621              matdiagprod)
622    formals(matdiagprod) <- alist(x=, y=NULL)
623    setMethod("tcrossprod", signature(x = "matrix", y = "diagonalMatrix"),
624              matdiagprod)
625    
626  setMethod("%*%", signature(x = "dgeMatrix", y = "diagonalMatrix"),  gediagprod <- function(x, y) {
           function(x, y) {  
627                dx <- dim(x)                dx <- dim(x)
628                dy <- dim(y)                dy <- dim(y)
629                if(dx[2] != dy[1]) stop("non-matching dimensions")                if(dx[2] != dy[1]) stop("non-matching dimensions")
630                if(y@diag == "N")                if(y@diag == "N")
631                    x@x <- x@x * rep(y@x, each = dx[1])                    x@x <- x@x * rep(y@x, each = dx[1])
632                x                x
633            })  }
634    setMethod("%*%", signature(x= "dgeMatrix", y= "diagonalMatrix"), gediagprod)
635    setMethod("%*%", signature(x= "lgeMatrix", y= "diagonalMatrix"), gediagprod)
636    formals(gediagprod) <- alist(x=, y=NULL)
637    setMethod("tcrossprod", signature(x = "dgeMatrix", y = "diagonalMatrix"),
638              gediagprod)
639    setMethod("tcrossprod", signature(x = "lgeMatrix", y = "diagonalMatrix"),
640              gediagprod)
641    
642  ## crossprod {more of these}  ## crossprod {more of these}
643    
644  ## tcrossprod --- all are not yet there: do the dense ones here:  ## tcrossprod --- all are not yet there: do the dense ones here:
645    
646    setMethod("%*%", signature(x = "diagonalMatrix", y = "denseMatrix"),
647              function(x, y) if(x@diag == "U") y else x %*% as(y, "generalMatrix"))
648    setMethod("%*%", signature(x = "denseMatrix", y = "diagonalMatrix"),
649              function(x, y) if(y@diag == "U") x else as(x, "generalMatrix") %*% y)
650    
651    
652  ## FIXME:  ## FIXME:
653  ## setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "denseMatrix"),  ## setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "denseMatrix"),
654  ##        function(x, y = NULL) {  ##        function(x, y = NULL) {
655  ##           })  ##           })
656    
657  ## setMethod("tcrossprod", signature(x = "denseMatrix", y = "diagonalMatrix"),  ##' @param x CsparseMatrix
658  ##        function(x, y = NULL) {  ##' @param y diagonalMatrix
659  ##           })  ##' @return x %*% y
660    Cspdiagprod <- function(x, y) {
661        dx <- dim(x <- .Call(Csparse_diagU2N, x))
662        dy <- dim(y)
663        if(dx[2] != dy[1]) stop("non-matching dimensions")
664        if(y@diag == "N") { ## otherwise: y == Diagonal(n) : multiplication is identity
665            if(!all(y@x[1L] == y@x[-1L]) && is(x, "symmetricMatrix"))
666                x <- as(x, "generalMatrix")
667            ind <- rep.int(seq_len(dx[2]), x@p[-1] - x@p[-dx[2]-1L])
668            x@x <- x@x * y@x[ind]
669            if(.hasSlot(x, "factors") && length(x@factors)) {# drop cashed ones
670                ## instead of dropping all factors, be smart about some
671                ## TODO ......
672                x@factors <- list()
673            }
674        }
675        x
676    }
677    
678    ##' @param x diagonalMatrix
679    ##' @param y CsparseMatrix
680    ##' @return x %*% y
681    diagCspprod <- function(x, y) {
682        dx <- dim(x)
683        dy <- dim(y <- .Call(Csparse_diagU2N, y))
684        if(dx[2] != dy[1]) stop("non-matching dimensions")
685        if(x@diag == "N") {
686            if(!all(x@x[1L] == x@x[-1L]) && is(y, "symmetricMatrix"))
687                y <- as(y, "generalMatrix")
688            y@x <- y@x * x@x[y@i + 1L]
689            if(.hasSlot(y, "factors") && length(yf <- y@factors)) {
690                ## TODO
691                if(FALSE) { ## instead of dropping all factors, be smart about some
692                    keep <- character()
693                    if(any(iLU <- names(yf) == "LU")) {
694                        keep <- "LU"
695                    }
696                    y@factors <- yf[keep]
697                } else y@factors <- list() ## for now
698            }
699        }
700        y
701    }
702    
703    setMethod("crossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
704              function(x, y = NULL) diagCspprod(x, y))
705    
706  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
707            function(x, y = NULL) { x <- as(x, "sparseMatrix"); callGeneric() })            function(x, y = NULL) diagCspprod(x, as(y, "CsparseMatrix")))
708    
709    ## Prefer calling diagCspprod to Cspdiagprod if going to transpose anyway
710    ##  x'y == (y'x)'
711    setMethod("crossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
712              function(x, y = NULL) t(diagCspprod(y, x)))
713    
714  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
715            function(x, y = NULL) { y <- as(y, "sparseMatrix"); callGeneric() })            function(x, y = NULL) t(diagCspprod(y, as(x, "Csparsematrix"))))
716    
717    setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
718              function(x, y = NULL) diagCspprod(x, t(y)))
719    
720  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
721            function(x, y = NULL) { x <- as(x, "sparseMatrix"); callGeneric() })            function(x, y = NULL) diagCspprod(x, t(as(y, "CsparseMatrix"))))
722    
723    setMethod("tcrossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
724              function(x, y = NULL) Cspdiagprod(x, y))
725    
726  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
727            function(x, y = NULL) { y <- as(y, "sparseMatrix"); callGeneric() })            function(x, y = NULL) Cspdiagprod(as(x, "CsparseMatrix"), y))
728    
729    setMethod("%*%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
730              function(x, y) diagCspprod(x, y))
731    
732    ## instead of "sparseMatrix", use: [RT]sparse.. ("closer" in method dispatch)
733    for(cl in c("TsparseMatrix", "RsparseMatrix")) {
734    
 ## FIXME?: In theory, this can be done *FASTER*, in some cases, via tapply1()  
735  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
736            function(x, y) as(x, "sparseMatrix") %*% y)            function(x, y) diagCspprod(as(x, "CsparseMatrix"), y))
 ## 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)  
 setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),  
           function(x, y) x %*% as(y, "CsparseMatrix"))  
 ## NB: this is *not* needed for Tsparse & Rsparse  
 ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*  
 ##       do indeed work by going through sparse (and *not* ddense)!  
737    
738  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
739            function(x, y) x %*% as(y, "sparseMatrix"))            function(x, y) Cspdiagprod(as(x, "CsparseMatrix"), y))
740    }
741    
742    setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
743              function(x, y) Cspdiagprod(x, y))
744    
745    ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*
746    ##       do indeed work by going through sparse (and *not* ddense)!
747    
748  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),
749            function(a, b, ...) {            function(a, b, ...) {
# Line 384  Line 753 
753            })            })
754    
755  solveDiag <- function(a, b, ...) {  solveDiag <- function(a, b, ...) {
756      if((n <- a@Dim[1]) != nrow(b))      if(a@Dim[1] != nrow(b))
757          stop("incompatible matrix dimensions")          stop("incompatible matrix dimensions")
758      ## trivially invert a 'in place' and multiply:      ## trivially invert a 'in place' and multiply:
759      a@x <- 1/ a@x      a@x <- 1/ a@x
# Line 396  Line 765 
765  setMethod("solve", signature(a = "diagonalMatrix", b = "Matrix"),  setMethod("solve", signature(a = "diagonalMatrix", b = "Matrix"),
766            solveDiag)            solveDiag)
767    
768    ## Schur()  ---> ./eigen.R
769    
770    
771    
772  ### ---------------- diagonal  o  sparse  -----------------------------  ###---------------- <Ops> (<Arith>, <Logic>, <Compare> ) ----------------------
   
773    
774  ## Use function for several signatures, in order to evade  ## Use function for several signatures, in order to evade
775  ## ambiguous dispatch for "ddi", since there's also Arith(ddense., ddense.)  diagOdiag <- function(e1,e2) {
776  diagOdiag <- function(e1,e2) { # result should also be diagonal      ## result should also be diagonal _ if possible _
777      r <- callGeneric(.diag.x(e1), .diag.x(e2)) # error if not "compatible"      r <- callGeneric(.diag.x(e1), .diag.x(e2)) # error if not "compatible"
778      if(is.numeric(r)) {      ## Check what happens with non-diagonals, i.e. (0 o 0), (FALSE o 0), ...:
779        r00 <- callGeneric(if(is.numeric(e1@x)) 0 else FALSE,
780                           if(is.numeric(e2@x)) 0 else FALSE)
781        if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* diagonal
782            if(is.numeric(r)) { # "double" *or* "integer"
783          if(is.numeric(e2@x)) {          if(is.numeric(e2@x)) {
784              e2@x <- r; return(.diag.2N(e2)) }              e2@x <- r; return(.diag.2N(e2)) }
785          if(!is.numeric(e1@x))          if(!is.numeric(e1@x))
786              ## e.g. e1, e2 are logical;              ## e.g. e1, e2 are logical;
787              e1 <- as(e1, "dMatrix")              e1 <- as(e1, "dMatrix")
788                if(!is.double(r)) r <- as.double(r)
789      }      }
790      else if(is.logical(r))      else if(is.logical(r))
791          e1 <- as(e1, "lMatrix")          e1 <- as(e1, "lMatrix")
792      else stop("intermediate 'r' is of type", typeof(r))          else stop(gettextf("intermediate 'r' is of type %s",
793                               typeof(r)), domain=NA)
794      e1@x <- r      e1@x <- r
795      .diag.2N(e1)      .diag.2N(e1)
796  }  }
797        else { ## result not diagonal, but at least symmetric:
798            ## e.g., m == m
799            isNum <- (is.numeric(r) || is.numeric(r00))
800            isLog <- (is.logical(r) || is.logical(r00))
801            Matrix.msg("exploding <diag> o <diag> into dense matrix", .M.level = 2)
802            d <- e1@Dim
803            n <- d[1]
804            stopifnot(length(r) == n)
805            if(isNum && !is.double(r)) r <- as.double(r)
806            xx <- as.vector(matrix(rbind(r, matrix(r00,n,n)), n,n))
807            newcl <-
808                paste0(if(isNum) "d" else if(isLog) {
809                    if(!any(is.na(r)) && !any(is.na(r00))) "n" else "l"
810                } else stop("not yet implemented .. please report"), "syMatrix")
811    
812            new(newcl, Dim = e1@Dim, Dimnames = e1@Dimnames, x = xx)
813        }
814    }
815    
816    ### This would be *the* way, but we get tons of "ambiguous method dispatch"
817    ## we use this hack instead of signature  x = "diagonalMatrix" :
818    diCls <- names(getClass("diagonalMatrix")@subclasses)
819    if(FALSE) {
820  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "diagonalMatrix"),  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "diagonalMatrix"),
821            diagOdiag)            diagOdiag)
822  ## These two are just for method disambiguation:  } else { ## These are just for method disambiguation:
823  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "diagonalMatrix"),      for(c1 in diCls)
824            diagOdiag)          for(c2 in diCls)
825  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "ddiMatrix"),              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)
826            diagOdiag)  }
827    
828    ## diagonal  o  triangular  |-->  triangular
829    ## diagonal  o  symmetric   |-->  symmetric
830    ##    {also when other is sparse: do these "here" --
831    ##     before conversion to sparse, since that loses "diagonality"}
832    diagOtri <- function(e1,e2) {
833        ## result must be triangular
834        r <- callGeneric(d1 <- .diag.x(e1), diag(e2)) # error if not "compatible"
835        ## Check what happens with non-diagonals, i.e. (0 o 0), (FALSE o 0), ...:
836        e1.0 <- if(.n1 <- is.numeric(d1   )) 0 else FALSE
837        r00 <- callGeneric(e1.0, if(.n2 <- is.numeric(e2[0])) 0 else FALSE)
838        if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* triangular
839            diag(e2) <- r
840            ## check what happens "in the triangle"
841            e2.2 <- if(.n2) 2 else TRUE
842            if(!callGeneric(e1.0, e2.2) == e2.2) { # values "in triangle" can change:
843                n <- dim(e2)[1L]
844                it <- indTri(n, upper = (e2@uplo == "U"))
845                e2[it] <- callGeneric(e1.0, e2[it])
846            }
847            e2
848        }
849        else { ## result not triangular ---> general
850            rr <- as(e2, "generalMatrix")
851            diag(rr) <- r
852            rr
853        }
854    }
855    
856    
857    setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "triangularMatrix"),
858              diagOtri)
859    ## For the reverse,  Ops == "Arith" | "Compare" | "Logic"
860    ##   'Arith'  :=  '"+"', '"-"', '"*"', '"^"', '"%%"', '"%/%"', '"/"'
861    setMethod("Arith", signature(e1 = "triangularMatrix", e2 = "diagonalMatrix"),
862              function(e1,e2)
863          { ## this must only trigger for *dense* e1
864              switch(.Generic,
865                     "+" = .Call(dtrMatrix_addDiag, as(e1,"dtrMatrix"),   .diag.x(e2)),
866                     "-" = .Call(dtrMatrix_addDiag, as(e1,"dtrMatrix"), - .diag.x(e2)),
867                     "*" = {
868                         n <- e2@Dim[1L]
869                         d2 <- if(e2@diag == "U") { # unit-diagonal
870                             d <- rep.int(as1(e2@x), n)
871                             e2@x <- d
872                             e2@diag <- "N"
873                             d
874                         } else e2@x
875                         e2@x <- diag(e1) * d2
876                         e2
877                     },
878                     "^" = { ## will be dense ( as  <ANY> ^ 0 == 1 ):
879                         e1 ^ as(e2, "denseMatrix")
880                     },
881                     ## otherwise:
882                     callGeneric(e1, diag2Tsmart(e2,e1)))
883    })
884    
885    ## Compare --> 'swap' (e.g.   e1 < e2   <==>  e2 > e1 ):
886    setMethod("Compare", signature(e1 = "triangularMatrix", e2 = "diagonalMatrix"),
887              .Cmp.swap)
888    ## '&' and "|'  are commutative:
889    setMethod("Logic", signature(e1 = "triangularMatrix", e2 = "diagonalMatrix"),
890              function(e1,e2) callGeneric(e2, e1))
891    
892  ## For almost everything else, diag* shall be treated "as sparse" :  ## For almost everything else, diag* shall be treated "as sparse" :
893  ## These are cheap implementations via coercion  ## These are cheap implementations via coercion
894    
895  ## for disambiguation  ## For disambiguation --- define this for "sparseMatrix" , then for "ANY";
896  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "sparseMatrix"),  ## and because we can save an .M.kind() call, we use this explicit
897            function(e1,e2) callGeneric(as(e1, "sparseMatrix"), e2))  ## "hack" for all diagonalMatrix *subclasses* instead of just "diagonalMatrix" :
898  setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "diagonalMatrix"),  ##
899            function(e1,e2) callGeneric(e1, as(e2, "sparseMatrix")))  ## ddi*:
900  ## in general:  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "sparseMatrix"),
901  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "ANY"),            function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2))
902            function(e1,e2) callGeneric(as(e1,"sparseMatrix"), e2))  setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "ddiMatrix"),
903  setMethod("Ops", signature(e1 = "ANY", e2 = "diagonalMatrix"),            function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "d")))
904            function(e1,e2) callGeneric(e1, as(e2,"sparseMatrix")))  ## ldi*
905    setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "sparseMatrix"),
906              function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2))
907    setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "ldiMatrix"),
908              function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "l")))
909    
910    ## Ops:  Arith  --> numeric : "dMatrix"
911    ##       Compare --> logical
912    ##       Logic   --> logical: "lMatrix"
913    
914    ## Other = "numeric" : stay diagonal if possible
915    ## ddi*: Arith: result numeric, potentially ddiMatrix
916    for(arg2 in c("numeric","logical"))
917    setMethod("Arith", signature(e1 = "ddiMatrix", e2 = arg2),
918              function(e1,e2) {
919                  n <- e1@Dim[1]
920                  f0 <- callGeneric(0, e2)
921                  if(all0(f0)) { # remain diagonal
922                      L1 <- (le <- length(e2)) == 1L
923                      if(e1@diag == "U") {
924                          if(any((r <- callGeneric(1, e2)) != 1)) {
925                              e1@diag <- "N"
926                              e1@x[seq_len(n)] <- r # possibly recycling r
927                          } ## else: result = e1  (is "U" diag)
928                      } else {
929                          r <- callGeneric(e1@x, e2)
930                          ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
931                          e1@x[] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
932                      }
933                      e1
934                  } else
935                      callGeneric(diag2tT.u(e1,e2, "d"), e2)
936              })
937    
938    for(arg1 in c("numeric","logical"))
939    setMethod("Arith", signature(e1 = arg1, e2 = "ddiMatrix"),
940              function(e1,e2) {
941                  n <- e2@Dim[1]
942                  f0 <- callGeneric(e1, 0)
943                  if(all0(f0)) { # remain diagonal
944                      L1 <- (le <- length(e1)) == 1L
945                      if(e2@diag == "U") {
946                          if(any((r <- callGeneric(e1, 1)) != 1)) {
947                              e2@diag <- "N"
948                              e2@x[seq_len(n)] <- r # possibly recycling r
949                          } ## else: result = e2  (is "U" diag)
950                      } else {
951                          r <- callGeneric(e1, e2@x)
952                          ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
953                          e2@x[] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
954                      }
955                      e2
956                  } else
957                      callGeneric(e1, diag2tT.u(e2,e1, "d"))
958              })
959    
960    ## ldi* Arith --> result numeric, potentially ddiMatrix
961    for(arg2 in c("numeric","logical"))
962    setMethod("Arith", signature(e1 = "ldiMatrix", e2 = arg2),
963              function(e1,e2) {
964                  n <- e1@Dim[1]
965                  f0 <- callGeneric(0, e2)
966                  if(all0(f0)) { # remain diagonal
967                      L1 <- (le <- length(e2)) == 1L
968                      E <- copyClass(e1, "ddiMatrix", c("diag", "Dim", "Dimnames"), check=FALSE)
969                      ## storage.mode(E@x) <- "double"
970                      if(e1@diag == "U") {
971                          if(any((r <- callGeneric(1, e2)) != 1)) {
972                              E@diag <- "N"
973                              E@x[seq_len(n)] <- r # possibly recycling r
974                          } ## else: result = E  (is "U" diag)
975                      } else {
976                          r <- callGeneric(e1@x, e2)
977                          ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
978                          E@x[seq_len(n)] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
979                      }
980                      E
981                  } else
982                      callGeneric(diag2tT.u(e1,e2, "l"), e2)
983              })
984    
985    for(arg1 in c("numeric","logical"))
986    setMethod("Arith", signature(e1 = arg1, e2 = "ldiMatrix"),
987              function(e1,e2) {
988                  n <- e2@Dim[1]
989                  f0 <- callGeneric(e1, 0)
990                  if(all0(f0)) { # remain diagonal
991                      L1 <- (le <- length(e1)) == 1L
992                      E <- copyClass(e2, "ddiMatrix", c("diag", "Dim", "Dimnames"), check=FALSE)
993                      ## storage.mode(E@x) <- "double"
994                      if(e2@diag == "U") {
995                          if(any((r <- callGeneric(e1, 1)) != 1)) {
996                              E@diag <- "N"
997                              E@x[seq_len(n)] <- r # possibly recycling r
998                          } ## else: result = E  (is "U" diag)
999                      } else {
1000                          r <- callGeneric(e1, e2@x)
1001                          ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
1002                          E@x[seq_len(n)] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
1003                      }
1004                      E
1005                  } else
1006                      callGeneric(e1, diag2tT.u(e2,e1, "l"))
1007              })
1008    
1009    ## ddi*: for "Ops" without "Arith": <Compare> or <Logic> --> result logical, potentially ldi
1010    ##
1011    ## Note that  ("numeric", "ddiMatrix")  is simply swapped, e.g.,
1012    if(FALSE) {
1013        selectMethod("<", c("numeric","lMatrix"))# Compare
1014        selectMethod("&", c("numeric","lMatrix"))# Logic
1015    } ## so we don't need to define a method here :
1016    
1017    for(arg2 in c("numeric","logical"))
1018    setMethod("Ops", signature(e1 = "ddiMatrix", e2 = arg2),
1019              function(e1,e2) {
1020                  n <- e1@Dim[1]
1021                  f0 <- callGeneric(0, e2)
1022                  if(all0(f0)) { # remain diagonal
1023                      L1 <- (le <- length(e2)) == 1L
1024                      E <- copyClass(e1, "ldiMatrix", c("diag", "Dim", "Dimnames"), check=FALSE)
1025                      ## storage.mode(E@x) <- "logical"
1026                      if(e1@diag == "U") {
1027                          if(any((r <- callGeneric(1, e2)) != 1)) {
1028                              E@diag <- "N"
1029                              E@x[seq_len(n)] <- r # possibly recycling r
1030                          } ## else: result = E  (is "U" diag)
1031                      } else {
1032                          r <- callGeneric(e1@x, e2)
1033                          ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
1034                          E@x[seq_len(n)] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
1035                      }
1036                      E
1037                  } else
1038                      callGeneric(diag2tT.u(e1,e2, "d"), e2)
1039              })
1040    
1041    ## ldi*: for "Ops" without "Arith": <Compare> or <Logic> --> result logical, potentially ldi
1042    for(arg2 in c("numeric","logical"))
1043    setMethod("Ops", signature(e1 = "ldiMatrix", e2 = arg2),
1044              function(e1,e2) {
1045                  n <- e1@Dim[1]
1046                  f0 <- callGeneric(FALSE, e2)
1047                  if(all0(f0)) { # remain diagonal
1048                      L1 <- (le <- length(e2)) == 1L
1049    
1050                      if(e1@diag == "U") {
1051                          if(any((r <- callGeneric(TRUE, e2)) != 1)) {
1052                              e1@diag <- "N"
1053                              e1@x[seq_len(n)] <- r # possibly recycling r
1054                          } ## else: result = e1  (is "U" diag)
1055                      } else {
1056                          r <- callGeneric(e1@x, e2)
1057                          ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
1058                          e1@x[] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
1059                      }
1060                      e1
1061                  } else
1062                      callGeneric(diag2tT.u(e1,e2, "l"), e2)
1063              })
1064    
1065    
1066    ## Not {"sparseMatrix", "numeric} :  {"denseMatrix", "matrix", ... }
1067    for(other in c("ANY", "Matrix", "dMatrix")) {
1068        ## ddi*:
1069        setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other),
1070                  function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2))
1071        setMethod("Ops", signature(e1 = other, e2 = "ddiMatrix"),
1072                  function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "d")))
1073        ## ldi*:
1074        setMethod("Ops", signature(e1 = "ldiMatrix", e2 = other),
1075                  function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2))
1076        setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),
1077                  function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "l")))
1078    }
1079    
1080    ## Direct subclasses of "denseMatrix": currently ddenseMatrix, ldense... :
1081    dense.subCl <- local({ dM.scl <- getClass("denseMatrix")@subclasses
1082                           names(dM.scl)[vapply(dM.scl, slot, 0, "distance") == 1] })
1083    for(DI in diCls) {
1084        dMeth <- if(extends(DI, "dMatrix"))
1085            function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2)
1086        else # "lMatrix", the only other kind for now
1087            function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2)
1088        for(c2 in c(dense.subCl, "Matrix")) {
1089            for(Fun in c("*", "&")) {
1090                setMethod(Fun, signature(e1 = DI, e2 = c2),
1091                          function(e1,e2) callGeneric(e1, Diagonal(x = diag(e2))))
1092                setMethod(Fun, signature(e1 = c2, e2 = DI),
1093                          function(e1,e2) callGeneric(Diagonal(x = diag(e1)), e2))
1094            }
1095            setMethod("^", signature(e1 = c2, e2 = DI),
1096                      function(e1,e2) callGeneric(Diagonal(x = diag(e1)), e2))
1097            for(Fun in c("%%", "%/%", "/")) ## 0 <op> 0 |--> NaN  for these.
1098                setMethod(Fun, signature(e1 = DI, e2 = c2), dMeth)
1099        }
1100    }
1101    
1102    
1103    ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"
1104    ### ----------   the last 4: separately here
1105    for(cl in diCls) {
1106    setMethod("any", cl,
1107              function (x, ..., na.rm) {
1108                  if(any(x@Dim == 0)) FALSE
1109                  else if(x@diag == "U") TRUE else any(x@x, ..., na.rm = na.rm)
1110              })
1111    setMethod("all",  cl, function (x, ..., na.rm) {
1112        n <- x@Dim[1]
1113        if(n >= 2) FALSE
1114        else if(n == 0 || x@diag == "U") TRUE
1115        else all(x@x, ..., na.rm = na.rm)
1116    })
1117    setMethod("prod", cl, function (x, ..., na.rm) {
1118        n <- x@Dim[1]
1119        if(n >= 2) 0
1120        else if(n == 0 || x@diag == "U") 1
1121        else ## n == 1, diag = "N" :
1122            prod(x@x, ..., na.rm = na.rm)
1123    })
1124    
1125    setMethod("sum", cl,
1126              function(x, ..., na.rm) {
1127                  r <- sum(x@x, ..., na.rm = na.rm)# double or integer, correctly
1128                  if(x@diag == "U" && !is.na(r)) r + x@Dim[1] else r
1129              })
1130    }
1131    
1132    ## The remaining ones are  max, min, range :
1133    
1134    setMethod("Summary", "ddiMatrix",
1135              function(x, ..., na.rm) {
1136                  if(any(x@Dim == 0)) callGeneric(numeric(0), ..., na.rm=na.rm)
1137                  else if(x@diag == "U")
1138                      callGeneric(x@x, 0, 1, ..., na.rm=na.rm)
1139                  else callGeneric(x@x, 0, ..., na.rm=na.rm)
1140              })
1141    setMethod("Summary", "ldiMatrix",
1142              function(x, ..., na.rm) {
1143                  if(any(x@Dim == 0)) callGeneric(logical(0), ..., na.rm=na.rm)
1144                  else if(x@diag == "U")
1145                      callGeneric(x@x, FALSE, TRUE, ..., na.rm=na.rm)
1146                  else callGeneric(x@x, FALSE, ..., na.rm=na.rm)
1147              })
1148    
1149    
1150    
# Line 455  Line 1159 
1159      invisible(x)      invisible(x)
1160  }  }
1161    
1162    ## somewhat consistent with "print" for sparseMatrix :
1163    setMethod("print", signature(x = "diagonalMatrix"), prDiag)
1164    
1165  setMethod("show", signature(object = "diagonalMatrix"),  setMethod("show", signature(object = "diagonalMatrix"),
1166            function(object) {            function(object) {
1167                d <- dim(object)                d <- dim(object)
1168                cl <- class(object)                cl <- class(object)
1169                cat(sprintf('%d x %d diagonal matrix of class "%s"\n',                cat(sprintf('%d x %d diagonal matrix of class "%s"',
1170                            d[1], d[2], cl))                            d[1], d[2], cl))
1171                  if(d[1] < 50) {
1172                      cat("\n")
1173                prDiag(object)                prDiag(object)
1174                  } else {
1175                      cat(", with diagonal entries\n")
1176                      show(diag(object))
1177                      invisible(object)
1178                  }
1179            })            })
1180    
1181    rm(arg1, arg2, other, DI, cl, c1, c2,
1182       dense.subCl, diCls)# not used elsewhere
1183    
1184    setMethod("summary", signature(object = "diagonalMatrix"),
1185              function(object, ...) {
1186                  d <- dim(object)
1187                  r <- summary(object@x, ...)
1188                  attr(r, "header") <-
1189                      sprintf('%d x %d diagonal Matrix of class "%s"',
1190                              d[1], d[2], class(object))
1191                  ## use ole' S3 technology for such a simple case
1192                  class(r) <- c("diagSummary", class(r))
1193                  r
1194              })
1195    
1196    print.diagSummary <- function (x, ...) {
1197        cat(attr(x, "header"),"\n")
1198        class(x) <- class(x)[-1]
1199        print(x, ...)
1200        invisible(x)
1201    }

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

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge