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 1708, Fri Dec 22 19:53:37 2006 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          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      ## slightly debatable if we really should return Csparse.. :  ## expand(<mer>) needed something like bdiag() for lower-triangular
117      as(ret, "CsparseMatrix")  ## (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  diag2T <- function(from) {          if(use.n) { ## return nsparseMatrix :
166      i <- if(from@diag == "U") integer(0) else seq_len(from@Dim[1]) - 1:1              r <- new("ntTMatrix")
167      new(paste(.M.kind(from), "tTMatrix", sep=''),          } 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)) {
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  ## is better than this:      n <- from@Dim[1]
215  ## setAs("diagonalMatrix", "sparseMatrix",      i <- seq_len(n) - 1L
216  ##       function(from)      new(paste0(kind, "sTMatrix"),
217  ##        as(from, if(is(from, "dMatrix")) "dgCMatrix" else "lgCMatrix"))          Dim = from@Dim, Dimnames = from@Dimnames,
218  setAs("diagonalMatrix", "CsparseMatrix",          i = i, j = i, uplo = uplo,
219        function(from) as(diag2T(from), "CsparseMatrix"))          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  setAs("diagonalMatrix", "matrix",  ## 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    diag2tT <- function(from) .diag2tT(from, "U", "d")
255    setAs("ddiMatrix", "triangularMatrix", diag2tT)
256    ##_no_longer_ setAs("ddiMatrix", "sparseMatrix", diag2tT)
257    ## needed too (otherwise <dense> -> Tsparse is taken):
258    setAs("ddiMatrix", "TsparseMatrix", diag2tT)
259    setAs("ddiMatrix", "dsparseMatrix", diag2tT)
260    setAs("ddiMatrix", "CsparseMatrix",
261          function(from) as(.diag2tT(from, "U", "d"), "CsparseMatrix"))
262    setAs("ddiMatrix", "symmetricMatrix",
263          function(from) .diag2sT(from, "U", "d"))
264    ##
265    ## ldi*:
266    diag2tT <- function(from) .diag2tT(from, "U", "l")
267    setAs("ldiMatrix", "triangularMatrix", diag2tT)
268    ##_no_longer_ setAs("ldiMatrix", "sparseMatrix", diag2tT)
269    ## needed too (otherwise <dense> -> Tsparse is taken):
270    setAs("ldiMatrix", "TsparseMatrix", diag2tT)
271    setAs("ldiMatrix", "lsparseMatrix", diag2tT)
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    
277    
278    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 163  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) {
379        x <- as(x, "CsparseMatrix") ## << was "TsparseMatrix" (Csparse is faster now)
380        x <- if(missing(i))
381            x[, j, drop=drop]
382        else if(missing(j))
383            if(nargs() == 4) x[i, , drop=drop] else x[i, drop=drop]
384        else
385            x[i,j, drop=drop]
386        if(isS4(x) && isDiagonal(x)) as(x, "diagonalMatrix") else x
387    }
388    
389    setMethod("[", signature(x = "diagonalMatrix", i = "index",
390                             j = "index", drop = "logical"), subDiag)
391    setMethod("[", signature(x = "diagonalMatrix", i = "index",
392                             j = "missing", drop = "logical"),
393              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",
401                             j = "index", drop = "logical"),
402              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    ## Only(?) current bug:  x[i] <- value  is wrong when  i is *vector*
408    replDiag <- function(x, i, j, ..., value) {
409        x <- as(x, "CsparseMatrix")# was "Tsparse.." till 2012-07
410        if(missing(i))
411            x[, j] <- value
412        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
417            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
423        if(isDiagonal(x)) as(x, "diagonalMatrix") else x
424    }
425    
426    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index",
427                                    j = "index", value = "replValue"), replDiag)
428    
429    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index",
430                                    j = "missing", value = "replValue"),
431                     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",
440                                    j = "index", value = "replValue"),
441                     function(x,i,j, ..., value) replDiag(x, j=j, value=value))
442    
443    ## x[] <- value :
444    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing",
445                                    j = "missing", value = "ANY"),
446                     function(x,i,j, ..., value)
447                 {
448                  if(all0(value)) { # be faster
449                      r <- new(paste0(.M.kindC(getClassDef(class(x))),"tTMatrix"))# of all "0"
450                      r@Dim <- x@Dim
451                      r@Dimnames <- x@Dimnames
452                      r
453                  } else { ## typically non-sense: assigning to full sparseMatrix
454                      x[TRUE] <- value
455                      x
456                  }
457              })
458    
459    
460  setReplaceMethod("[", signature(x = "diagonalMatrix",  setReplaceMethod("[", signature(x = "diagonalMatrix",
461                                  i = "ANY", j = "ANY", value = "ANY"),                                  i = "matrix", # 2-col.matrix
462                   function(x, i, j, value) {                                  j = "missing", value = "replValue"),
463                       r <- callGeneric(x = as(x,"sparseMatrix"),                   function(x,i,j, ..., value) {
464                                        i=i, j=j, value=value)                       if(ncol(i) == 2) {
465                       if(isDiagonal(r)) as(r, "diagonalMatrix") else r                           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      }      }
     x  
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 250  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 273  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"),  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  ### ---------------- diagonal  o  sparse  -----------------------------  setMethod("crossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
697              function(x, y = NULL) diagCspprod(x, y))
698    
699    setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
700              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"),
708              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"),
714              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"),
720              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    
725    setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
726              function(x, y) diagCspprod(as(x, "CsparseMatrix"), y))
727    
728    setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
729              function(x, y) Cspdiagprod(as(x, "CsparseMatrix"), y))
730    
731    setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
732              function(x, y) Cspdiagprod(x, y))
733    
734    ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*
735    ##       do indeed work by going through sparse (and *not* ddense)!
736    
737    setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),
738              function(a, b, ...) {
739                  a@x <- 1/ a@x
740                  a@Dimnames <- a@Dimnames[2:1]
741                  a
742              })
743    
744    solveDiag <- function(a, b, ...) {
745        if(a@Dim[1] != nrow(b))
746            stop("incompatible matrix dimensions")
747        ## trivially invert a 'in place' and multiply:
748        a@x <- 1/ a@x
749        a@Dimnames <- a@Dimnames[2:1]
750        a %*% b
751    }
752    setMethod("solve", signature(a = "diagonalMatrix", b = "matrix"),
753              solveDiag)
754    setMethod("solve", signature(a = "diagonalMatrix", b = "Matrix"),
755              solveDiag)
756    
757    ## Schur()  ---> ./eigen.R
758    
759    
760    
761    ###---------------- <Ops> (<Arith>, <Logic>, <Compare> ) ----------------------
762    
763  ## Use function for several signatures, in order to evade  ## Use function for several signatures, in order to evade
764  ## ambiguous dispatch for "ddi", since there's also Arith(ddense., ddense.)  diagOdiag <- function(e1,e2) {
765  diagOdiag <- function(e1,e2) { # result should also be diagonal      ## 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"
767      if(is.numeric(r)) {      ## Check what happens with non-diagonals, i.e. (0 o 0), (FALSE o 0), ...:
768        r00 <- callGeneric(if(is.numeric(e1@x)) 0 else FALSE,
769                           if(is.numeric(e2@x)) 0 else FALSE)
770        if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* diagonal
771            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:
787            ## e.g., m == m
788            isNum <- (is.numeric(r) || is.numeric(r00))
789            isLog <- (is.logical(r) || is.logical(r00))
790            Matrix.msg("exploding <diag> o <diag> into dense matrix", .M.level = 2)
791            d <- e1@Dim
792            n <- d[1]
793            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))
796            newcl <-
797                paste0(if(isNum) "d" else if(isLog) {
798                    if(!any(is.na(r)) && !any(is.na(r00))) "n" else "l"
799                } else stop("not yet implemented .. please report"), "syMatrix")
800    
801            new(newcl, Dim = e1@Dim, Dimnames = e1@Dimnames, x = xx)
802        }
803    }
804    
805    ### This would be *the* way, but we get tons of "ambiguous method dispatch"
806    ## we use this hack instead of signature  x = "diagonalMatrix" :
807    diCls <- names(getClass("diagonalMatrix")@subclasses)
808    if(FALSE) {
809  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "diagonalMatrix"),  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "diagonalMatrix"),
810            diagOdiag)            diagOdiag)
811  ## These two are just for method disambiguation:  } else { ## These are just for method disambiguation:
812  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "diagonalMatrix"),      for(c1 in diCls)
813            diagOdiag)          for(c2 in diCls)
814  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "ddiMatrix"),              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)
815            diagOdiag)  }
816    
817    ## diagonal  o  triangular  |-->  triangular
818    ## diagonal  o  symmetric   |-->  symmetric
819    ##    {also when other is sparse: do these "here" --
820    ##     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
883    
884  ## for disambiguation  ## For disambiguation --- define this for "sparseMatrix" , then for "ANY";
885  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "sparseMatrix"),  ## and because we can save an .M.kind() call, we use this explicit
886            function(e1,e2) callGeneric(as(e1, "sparseMatrix"), e2))  ## "hack" for all diagonalMatrix *subclasses* instead of just "diagonalMatrix" :
887  setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "diagonalMatrix"),  ##
888            function(e1,e2) callGeneric(e1, as(e2, "sparseMatrix")))  ## ddi*:
889  ## in general:  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "sparseMatrix"),
890  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "ANY"),            function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2))
891            function(e1,e2) callGeneric(as(e1,"sparseMatrix"), e2))  setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "ddiMatrix"),
892  setMethod("Ops", signature(e1 = "ANY", e2 = "diagonalMatrix"),            function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "d")))
893            function(e1,e2) callGeneric(e1, as(e2,"sparseMatrix")))  ## ldi*
894    setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "sparseMatrix"),
895              function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2))
896    setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "ldiMatrix"),
897              function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "l")))
898    
899    ## Ops:  Arith  --> numeric : "dMatrix"
900    ##       Compare --> logical
901    ##       Logic   --> logical: "lMatrix"
902    
903    ## Other = "numeric" : stay diagonal if possible
904    ## ddi*: Arith: result numeric, potentially ddiMatrix
905    for(arg2 in c("numeric","logical"))
906    setMethod("Arith", signature(e1 = "ddiMatrix", e2 = arg2),
907              function(e1,e2) {
908                  n <- e1@Dim[1]
909                  f0 <- callGeneric(0, e2)
910                  if(all(is0(f0))) { # remain diagonal
911                      L1 <- (le <- length(e2)) == 1L
912                      if(e1@diag == "U") {
913                          if(any((r <- callGeneric(1, e2)) != 1)) {
914                              e1@diag <- "N"
915                              e1@x[seq_len(n)] <- r # possibly recycling r
916                          } ## else: result = e1  (is "U" diag)
917                      } else {
918                          r <- callGeneric(e1@x, e2)
919                          ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
920                          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)
925              })
926    
927    for(arg1 in c("numeric","logical"))
928    setMethod("Arith", signature(e1 = arg1, e2 = "ddiMatrix"),
929              function(e1,e2) {
930                  n <- e2@Dim[1]
931                  f0 <- callGeneric(e1, 0)
932                  if(all(is0(f0))) { # remain diagonal
933                      L1 <- (le <- length(e1)) == 1L
934                      if(e2@diag == "U") {
935                          if(any((r <- callGeneric(e1, 1)) != 1)) {
936                              e2@diag <- "N"
937                              e2@x[seq_len(n)] <- r # possibly recycling r
938                          } ## else: result = e2  (is "U" diag)
939                      } else {
940                          r <- callGeneric(e1, e2@x)
941                          ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
942                          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"))
947              })
948    
949    ## ldi* Arith --> result numeric, potentially ddiMatrix
950    for(arg2 in c("numeric","logical"))
951    setMethod("Arith", signature(e1 = "ldiMatrix", e2 = arg2),
952              function(e1,e2) {
953                  n <- e1@Dim[1]
954                  f0 <- callGeneric(0, e2)
955                  if(all(is0(f0))) { # remain diagonal
956                      L1 <- (le <- length(e2)) == 1L
957                      E <- copyClass(e1, "ddiMatrix", c("diag", "Dim", "Dimnames"))#FIXME: if ok, check=FALSE
958                      ## E <- copyClass(e1, "ddiMatrix", check=FALSE)
959                      ## storage.mode(E@x) <- "double"
960                      if(e1@diag == "U") {
961                          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)
967                          ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
968                          E@x[seq_len(n)] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
969                      }
970                      E
971                  } else
972                      callGeneric(diag2tT.u(e1,e2, "l"), e2)
973              })
974    
975  ## FIXME?: In theory, this can be done *FASTER*, in some cases, via tapply1()  for(arg1 in c("numeric","logical"))
976  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("Arith", signature(e1 = arg1, e2 = "ldiMatrix"),
977            function(x, y) as(x, "sparseMatrix") %*% y)            function(e1,e2) {
978  ## NB: The previous is *not* triggering for  "ddi" o "dgC" (= distance 3)                n <- e2@Dim[1]
979  ##     since there's a "ddense" o "Csparse" at dist. 2 => triggers first.                f0 <- callGeneric(e1, 0)
980  ## ==> do this:                if(all(is0(f0))) { # remain diagonal
981  setMethod("%*%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),                    L1 <- (le <- length(e1)) == 1L
982            function(x, y) as(x, "CsparseMatrix") %*% y)                    E <- copyClass(e2, "ddiMatrix", c("diag", "Dim", "Dimnames"))#FIXME: if ok, check=FALSE
983  ## NB: this is *not* needed for Tsparse & Rsparse                    ## E <- copyClass(e2, "ddiMatrix", check=FALSE)
984  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*                    ## storage.mode(E@x) <- "double"
985  ##       do indeed work by going throug sparse (and *not* ddense)!                    if(e2@diag == "U") {
986                          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)
992                          ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
993                          E@x[seq_len(n)] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
994                      }
995                      E
996                  } else
997                      callGeneric(e1, diag2tT.u(e2,e1, "l"))
998              })
999    
1000  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),  ## ddi*: for "Ops" without "Arith": <Compare> or <Logic> --> result logical, potentially ldi
1001            function(x, y) x %*% as(y, "sparseMatrix"))  ##
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) {
1011                  n <- e1@Dim[1]
1012                  f0 <- callGeneric(0, e2)
1013                  if(all(is0(f0))) { # remain diagonal
1014                      L1 <- (le <- length(e2)) == 1L
1015                      E <- copyClass(e1, "ldiMatrix", c("diag", "Dim", "Dimnames"))#FIXME: if ok, check=FALSE
1016                      ## E <- copyClass(e1, "ldiMatrix", check=FALSE)
1017                      ## storage.mode(E@x) <- "logical"
1018                      if(e1@diag == "U") {
1019                          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)
1025                          ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
1026                          E@x[seq_len(n)] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
1027                      }
1028                      E
1029                  } else
1030                      callGeneric(diag2tT.u(e1,e2, "d"), e2)
1031              })
1032    
1033  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  ## ldi*: for "Ops" without "Arith": <Compare> or <Logic> --> result logical, potentially ldi
1034            function(x, y = NULL) { x <- as(x, "sparseMatrix"); callGeneric() })  for(arg2 in c("numeric","logical"))
1035    setMethod("Ops", signature(e1 = "ldiMatrix", e2 = arg2),
1036              function(e1,e2) {
1037                  n <- e1@Dim[1]
1038                  f0 <- callGeneric(FALSE, e2)
1039                  if(all(is0(f0))) { # remain diagonal
1040                      L1 <- (le <- length(e2)) == 1L
1041    
1042  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),                    if(e1@diag == "U") {
1043            function(x, y = NULL) { y <- as(y, "sparseMatrix"); callGeneric() })                        if(any((r <- callGeneric(TRUE, e2)) != 1)) {
1044                              e1@diag <- "N"
1045                              e1@x[seq_len(n)] <- r # possibly recycling r
1046                          } ## else: result = e1  (is "U" diag)
1047                      } else {
1048                          r <- callGeneric(e1@x, e2)
1049                          ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
1050                          e1@x[] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
1051                      }
1052                      e1
1053                  } else
1054                      callGeneric(diag2tT.u(e1,e2, "l"), e2)
1055              })
1056    
 setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  
           function(x, y = NULL) { x <- as(x, "sparseMatrix"); callGeneric() })  
1057    
1058  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  ## Not {"sparseMatrix", "numeric} :  {"denseMatrix", "matrix", ... }
1059            function(x, y = NULL) { y <- as(y, "sparseMatrix"); callGeneric() })  for(other in c("ANY", "Matrix", "dMatrix")) {
1060        ## ddi*:
1061        setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other),
1062                  function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2))
1063        setMethod("Ops", signature(e1 = other, e2 = "ddiMatrix"),
1064                  function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "d")))
1065        ## ldi*:
1066        setMethod("Ops", signature(e1 = "ldiMatrix", e2 = other),
1067                  function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2))
1068        setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),
1069                  function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "l")))
1070    }
1071    
1072    ## Direct subclasses of "denseMatrix": currently ddenseMatrix, ldense... :
1073    dense.subCl <- local({ dM.scl <- getClass("denseMatrix")@subclasses
1074                           names(dM.scl)[vapply(dM.scl, slot, 0, "distance") == 1] })
1075    for(DI in diCls) {
1076        dMeth <- if(extends(DI, "dMatrix"))
1077            function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2)
1078        else # "lMatrix", the only other kind for now
1079            function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2)
1080        for(c2 in c(dense.subCl, "Matrix")) {
1081            for(Fun in c("*", "&")) {
1082                setMethod(Fun, signature(e1 = DI, e2 = c2),
1083                          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"
1096    ### ----------   the last 4: separately here
1097    for(cl in diCls) {
1098    setMethod("any", cl,
1099              function (x, ..., na.rm) {
1100                  if(any(x@Dim == 0)) FALSE
1101                  else if(x@diag == "U") TRUE else any(x@x, ..., na.rm = na.rm)
1102              })
1103    setMethod("all",  cl, function (x, ..., na.rm) {
1104        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,
1118              function(x, ..., na.rm) {
1119                  r <- sum(x@x, ..., na.rm = na.rm)# double or integer, correctly
1120                  if(x@diag == "U" && !is.na(r)) r + x@Dim[1] else r
1121              })
1122    }
1123    
1124    ## The remaining ones are  max, min, range :
1125    
1126    setMethod("Summary", "ddiMatrix",
1127              function(x, ..., na.rm) {
1128                  if(any(x@Dim == 0)) callGeneric(numeric(0), ..., na.rm=na.rm)
1129                  else if(x@diag == "U")
1130                      callGeneric(x@x, 0, 1, ..., na.rm=na.rm)
1131                  else callGeneric(x@x, 0, ..., na.rm=na.rm)
1132              })
1133    setMethod("Summary", "ldiMatrix",
1134              function(x, ..., na.rm) {
1135                  if(any(x@Dim == 0)) callGeneric(logical(0), ..., na.rm=na.rm)
1136                  else if(x@diag == "U")
1137                      callGeneric(x@x, FALSE, TRUE, ..., na.rm=na.rm)
1138                  else callGeneric(x@x, FALSE, ..., na.rm=na.rm)
1139              })
1140    
1141    
1142    
1143  ## similar to prTriang() in ./Auxiliaries.R :  ## similar to prTriang() in ./Auxiliaries.R :
# Line 395  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.1708  
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