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 2840, Fri Oct 5 22:18:37 2012 UTC
# Line 5  Line 5 
5  ##          but *not* diag() extractor!  ##          but *not* diag() extractor!
6  Diagonal <- function(n, x = NULL)  Diagonal <- function(n, x = NULL)
7  {  {
8      ## Allow  Diagonal(4)  and  Diagonal(x=1:5)      ## Allow  Diagonal(4), Diagonal(x=1:5), and  Diagonal(4, TRUE)
9      if(missing(n))      if(missing(n))
10          n <- length(x)          n <- length(x)
11      else {      else {
# Line 16  Line 16 
16      if(missing(x)) ## unit diagonal matrix      if(missing(x)) ## unit diagonal matrix
17          new("ddiMatrix", Dim = c(n,n), diag = "U")          new("ddiMatrix", Dim = c(n,n), diag = "U")
18      else {      else {
19          stopifnot(length(x) == n)          lx <- length(x)
20            lx.1 <- lx == 1L
21            stopifnot(lx.1 || lx == n) # but keep 'x' short for now
22          if(is.logical(x))          if(is.logical(x))
23              cl <- "ldiMatrix"              cl <- "ldiMatrix"
24          else if(is.numeric(x)) {          else if(is.numeric(x)) {
# Line 26  Line 28 
28          else if(is.complex(x)) {          else if(is.complex(x)) {
29              cl <- "zdiMatrix"  # will not yet work              cl <- "zdiMatrix"  # will not yet work
30          } else stop("'x' has invalid data type")          } else stop("'x' has invalid data type")
31          new(cl, Dim = c(n,n), diag = "N", x = x)          if(lx.1 && !is.na(x) && x == 1) # cheap check for uni-diagonal..
32                new(cl, Dim = c(n,n), diag = "U")
33            else
34                new(cl, Dim = c(n,n), diag = "N",
35                    x = if(lx.1) rep.int(x,n) else x)
36      }      }
37  }  }
38    
39  ### This is modified from a post of Bert Gunter to R-help on  1 Sep 2005.  .sparseDiagonal <- function(n, x = 1, uplo = "U",
40  ### Bert's code built on a post by Andy Liaw who most probably was influenced                              shape = if(missing(cols)) "t" else "g",
41  ### by earlier posts, notably one by Scott Chasalow on S-news, 16 Jan 2002                              unitri, kind,
42  ### who posted his bdiag() function written in December 1995.                              cols = if(n) 0:(n - 1L) else integer(0))
43    {
44        stopifnot(n == (n. <- as.integer(n)), (n <- n.) >= 0)
45        if(!(mcols <- missing(cols)))
46            stopifnot(0 <= (cols <- as.integer(cols)), cols < n)
47        m <- length(cols)
48        if(missing(kind))
49            kind <-
50                if(is.double(x)) "d"
51                else if(is.logical(x)) "l"
52                else { ## for now
53                    storage.mode(x) <- "double"
54                    "d"
55                }
56        else stopifnot(any(kind == c("d","l","n")))
57        stopifnot(is.character(shape), nchar(shape) == 1,
58                  any(shape == c("t","s","g"))) # triangular / symmetric / general
59        if((missing(unitri) || unitri) && shape == "t" &&
60           (mcols || cols == 0:(n-1L)) &&
61           ((any(kind == c("l", "n")) && allTrue(x)) ||
62            (    kind == "d"          && allTrue(x == 1)))) { ## uni-triangular
63            new(paste0(kind,"tCMatrix"), Dim = c(n,n),
64                       uplo = uplo, diag = "U", p = rep.int(0L, n+1L))
65        }
66        else if(kind == "n") {
67            if(shape == "g")
68                new("ngCMatrix", Dim = c(n,m), i = cols, p = 0:m)
69            else new(paste0("n", shape, "CMatrix"), Dim = c(n,m), uplo = uplo,
70                     i = cols, p = 0:m)
71        }
72        else { ## kind != "n" -- have x slot :
73            if((lx <- length(x)) == 1) x <- rep.int(x, m)
74            else if(lx != m) stop("length(x) must be either 1 or #{cols}")
75            if(shape == "g")
76                new(paste0(kind, "gCMatrix"), Dim = c(n,m),
77                    x = x, i = cols, p = 0:m)
78            else new(paste0(kind, shape, "CMatrix"), Dim = c(n,m), uplo = uplo,
79                     x = x, i = cols, p = 0:m)
80        }
81    }
82    
83  bdiag <- function(...) {  ## Pkg 'spdep' had (relatively slow) versions of this as_dsCMatrix_I()
84      if(nargs() == 0) return(new("dgCMatrix"))  .symDiagonal <- function(n, x = rep.int(1,n), uplo = "U")
85      ## else :      .sparseDiagonal(n, x, uplo, shape = "s")
86      mlist <- if (nargs() == 1) as.list(...) else list(...)  
87      dims <- sapply(mlist, dim)  # instead of   diagU2N(as(Diagonal(n), "CsparseMatrix")), diag = "N" in any case:
88    .trDiagonal <- function(n, x = 1, uplo = "U", unitri=TRUE)
89        .sparseDiagonal(n, x, uplo, shape = "t", unitri=unitri)
90    
91    
92    ## This is modified from a post of Bert Gunter to R-help on  1 Sep 2005.
93    ## Bert's code built on a post by Andy Liaw who most probably was influenced
94    ## by earlier posts, notably one by Scott Chasalow on S-news, 16 Jan 2002
95    ## who posted his bdiag() function written in December 1995.
96    if(FALSE)##--- no longer used:
97    .bdiag <- function(lst) {
98        ## block-diagonal matrix [a dgTMatrix] from list of matrices
99        stopifnot(is.list(lst), length(lst) >= 1)
100        dims <- vapply(lst, dim, 1L, USE.NAMES=FALSE)
101      ## make sure we had all matrices:      ## make sure we had all matrices:
102      if(!(is.matrix(dims) && nrow(dims) == 2))      if(!(is.matrix(dims) && nrow(dims) == 2))
103          stop("some arguments are not matrices")          stop("some arguments are not matrices")
104      csdim <- rbind(rep.int(0:0, 2),      csdim <- rbind(rep.int(0L, 2),
105                     apply(sapply(mlist, dim), 1, cumsum))                     apply(dims, 1, cumsum))
106      ret <- new("dgTMatrix", Dim = as.integer(csdim[nrow(csdim),]))      r <- new("dgTMatrix")
107        r@Dim <- as.integer(csdim[nrow(csdim),])
108      add1 <- matrix(1:0, 2,2)      add1 <- matrix(1:0, 2,2)
109      for(i in seq_along(mlist)) {      for(i in seq_along(lst)) {
110          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])
111          if(is.null(dim(indx))) ## non-square matrix          if(is.null(dim(indx))) ## non-square matrix
112              ret[indx[[1]],indx[[2]]] <- mlist[[i]]              r[indx[[1]],indx[[2]]] <- lst[[i]]
113          else ## square matrix          else ## square matrix
114              ret[indx[,1],indx[,2]] <- mlist[[i]]              r[indx[,1], indx[,2]] <- lst[[i]]
115        }
116        r
117    }
118    ## expand(<mer>) needed something like bdiag() for lower-triangular
119    ## (Tsparse) Matrices; hence Doug Bates provided a much more efficient
120    ##  implementation for those; now extended and generalized:
121    .bdiag <- function(lst) {
122        ## block-diagonal matrix [a dgTMatrix] from list of matrices
123        stopifnot(is.list(lst), (nl <- length(lst)) >= 1)
124    
125        Tlst <- lapply(lapply(lst, Matrix:::as_Csp2), # includes "diagU2N"
126                       as, "TsparseMatrix")
127        if(nl == 1) return(Tlst[[1]])
128        ## else
129        i_off <- c(0L, cumsum(vapply(Tlst, nrow, 1L)))
130        j_off <- c(0L, cumsum(vapply(Tlst, ncol, 1L)))
131    
132        clss <- vapply(Tlst, class, "")
133        typ <- substr(clss, 2, 2)
134        knd <- substr(clss, 1, 1)
135        sym <- typ == "s" # symmetric ones
136        tri <- typ == "t" # triangular ones
137        use.n <- any(is.n <- knd == "n")
138        if(use.n && !(use.n <- all(is.n))) {
139            Tlst[is.n] <- lapply(Tlst[is.n], as, "lMatrix")
140            knd [is.n] <- "l"
141        }
142        use.l <- !use.n && all(knd == "l")
143        if(all(sym)) { ## result should be *symmetric*
144            uplos <- vapply(Tlst, slot, ".", "uplo") ## either "U" or "L"
145            tLU <- table(uplos)# of length 1 or 2 ..
146            if(length(tLU) == 1) { ## all "U" or all "L"
147                useU <- uplos[1] == "U"
148            } else { ## length(tLU) == 2, counting "L" and "U"
149                useU <- diff(tLU) >= 0
150                if(useU && (hasL <- tLU[1] > 0))
151                    Tlst[hasL] <- lapply(Tlst[hasL], t)
152                else if(!useU && (hasU <- tLU[2] > 0))
153                    Tlst[hasU] <- lapply(Tlst[hasU], t)
154            }
155            if(use.n) { ## return nsparseMatrix :
156                r <- new("nsTMatrix")
157            } else {
158                r <- new(paste0(if(use.l) "l" else "d", "sTMatrix"))
159                r@x <- unlist(lapply(Tlst, slot, "x"))
160      }      }
161      ## slightly debatable if we really should return Csparse.. :          r@uplo <- if(useU) "U" else "L"
     as(ret, "CsparseMatrix")  
162  }  }
163        else if(all(tri) && { ULs <- vapply(Tlst, slot, ".", "uplo")##  "U" or "L"
164                              all(ULs[1L] == ULs[-1L]) } ## all upper or all lower
165           ){ ## *triangular* result
166    
167  diag2T <- function(from) {          if(use.n) { ## return nsparseMatrix :
168      i <- if(from@diag == "U") integer(0) else seq_len(from@Dim[1]) - 1:1              r <- new("ntTMatrix")
169      new(paste(.M.kind(from), "tTMatrix", sep=''),          } else {
170                r <- new(paste0(if(use.l) "l" else "d", "tTMatrix"))
171                r@x <- unlist(lapply(Tlst, slot, "x"))
172            }
173            r@uplo <- ULs[1L]
174        }
175        else {
176            if(any(sym))
177                Tlst[sym] <- lapply(Tlst[sym], as, "generalMatrix")
178            if(use.n) { ## return nsparseMatrix :
179                r <- new("ngTMatrix")
180            } else {
181                r <- new(paste0(if(use.l) "l" else "d", "gTMatrix"))
182                r@x <- unlist(lapply(Tlst, slot, "x"))
183            }
184        }
185        r@Dim <- c(i_off[nl+1], j_off[nl + 1])
186        r@i <- unlist(lapply(1:nl, function(k) Tlst[[k]]@i + i_off[k]))
187        r@j <- unlist(lapply(1:nl, function(k) Tlst[[k]]@j + j_off[k]))
188        r
189    }
190    
191    bdiag <- function(...) {
192        if((nA <- nargs()) == 0) return(new("dgCMatrix"))
193        if(nA == 1 && !is.list(...))
194            return(as(..., "CsparseMatrix"))
195        alis <- if(nA == 1 && is.list(..1)) ..1 else list(...)
196        if(length(alis) == 1)
197            return(as(alis[[1]], "CsparseMatrix"))
198    
199        ## else : two or more arguments
200        as(.bdiag(alis), "CsparseMatrix")
201    }
202    
203    
204    .diag2tT <- function(from, uplo = "U", kind = .M.kind(from)) {
205        ## to triangular Tsparse
206        i <- if(from@diag == "U") integer(0) else seq_len(from@Dim[1]) - 1L
207        new(paste0(kind, "tTMatrix"),
208          diag = from@diag, Dim = from@Dim, Dimnames = from@Dimnames,          diag = from@diag, Dim = from@Dim, Dimnames = from@Dimnames,
209            uplo = uplo,
210          x = from@x, # <- ok for diag = "U" and "N" (!)          x = from@x, # <- ok for diag = "U" and "N" (!)
211          i = i, j = i)          i = i, j = i)
212  }  }
213    
214  setAs("diagonalMatrix", "triangularMatrix", diag2T)  .diag2sT <- function(from, uplo = "U", kind = .M.kind(from)) {
215  setAs("diagonalMatrix", "sparseMatrix", diag2T)      ## to symmetric Tsparse
216        n <- from@Dim[1]
217        i <- seq_len(n) - 1L
218        new(paste0(kind, "sTMatrix"),
219            Dim = from@Dim, Dimnames = from@Dimnames,
220            i = i, j = i, uplo = uplo,
221            x = if(from@diag == "N") from@x else ## "U"-diag
222            rep.int(switch(kind,
223                           "d" = 1.,
224                           "l" =,
225                           "n" = TRUE,
226                           ## otherwise
227                           stop("'", kind,"' kind not yet implemented")), 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    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):  ## needed too (otherwise <dense> -> Tsparse is taken):
258  setAs("diagonalMatrix", "TsparseMatrix", diag2T)  setAs("ddiMatrix", "TsparseMatrix", diag2tT)
259  ## is better than this:  setAs("ddiMatrix", "dsparseMatrix", diag2tT)
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    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  setAs("diagonalMatrix", "matrix",  
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) if(m@diag == "U") rep.int(as1(m@x), m@Dim[1]) else m@x
320    
321    .diag.2N <- function(m) {
322        if(m@diag == "U") m@diag <- "N"
323        m
324    }
325    
 if(FALSE) # now have faster  "ddense" -> "dge"  
326  setAs("ddiMatrix", "dgeMatrix",  setAs("ddiMatrix", "dgeMatrix",
327        function(from) as(as(from, "matrix"), "dgeMatrix"))        function(from) .Call(dup_mMatrix_as_dgeMatrix, from))
328    setAs("ddiMatrix", "ddenseMatrix",
329          function(from) as(as(from, "triangularMatrix"),"denseMatrix"))
330    setAs("ldiMatrix", "ldenseMatrix",
331          function(from) as(as(from, "triangularMatrix"),"denseMatrix"))
332    
333    
334  setAs("matrix", "diagonalMatrix",  setAs("matrix", "diagonalMatrix",
335        function(from) {        function(from) {
336            d <- dim(from)            d <- dim(from)
337            if(d[1] != (n <- d[2])) stop("non-square matrix")            if(d[1] != (n <- d[2])) stop("non-square matrix")
338            if(any(from[row(from) != col(from)] != 0))            if(any(from[row(from) != col(from)] != 0))
339                stop("has non-zero off-diagonal entries")                stop("matrix with non-zero off-diagonals cannot be coerced to diagonalMatrix")
340            x <- diag(from)            x <- diag(from)
341            if(is.logical(x)) {            if(is.logical(x)) {
342                cl <- "ldiMatrix"                cl <- "ldiMatrix"
343                uni <- all(x)                uni <- allTrue(x) ## uni := {is it unit-diagonal ?}
344            } else {            } else {
345                cl <- "ddiMatrix"                cl <- "ddiMatrix"
346                uni <- all(x == 1)                uni <- allTrue(x == 1)
347                storage.mode(x) <- "double"                storage.mode(x) <- "double"
348            } ## TODO: complex            } ## TODO: complex
349            new(cl, Dim = c(n,n), diag = if(uni) "U" else "N",            new(cl, Dim = c(n,n), diag = if(uni) "U" else "N",
# Line 165  Line 360 
360            x <- diag(from)            x <- diag(from)
361            if(is.logical(x)) {            if(is.logical(x)) {
362                cl <- "ldiMatrix"                cl <- "ldiMatrix"
363                uni <- all(x)                uni <- allTrue(x)
364            } else {            } else {
365                cl <- "ddiMatrix"                cl <- "ddiMatrix"
366                uni <- all(x == 1)                uni <- allTrue(x == 1)
367                storage.mode(x) <- "double"                storage.mode(x) <- "double"
368            }            } ## TODO: complex
369            new(cl, Dim = c(n,n), diag = if(uni) "U" else "N",            new(cl, Dim = c(n,n), diag = if(uni) "U" else "N",
370                x = if(uni) x[FALSE] else x)                x = if(uni) x[FALSE] else x)
371        })        })
372    
373    
374  setMethod("diag", signature(x = "diagonalMatrix"),  setMethod("diag", signature(x = "diagonalMatrix"),
375            function(x = 1, nrow, ncol = n) .diag.x(x))            function(x = 1, nrow, ncol) .diag.x(x))
   
376    
377  subDiag <- function(x, i, j, drop) {  subDiag <- function(x, i, j, ..., drop) {
378      x <- as(x, "sparseMatrix")      x <- as(x, "CsparseMatrix") ## << was "TsparseMatrix" (Csparse is faster now)
379      x <- if(missing(i))      x <- if(missing(i))
380          x[, j, drop=drop]          x[, j, drop=drop]
381      else if(missing(j))      else if(missing(j))
382          x[i, , drop=drop]          if(nargs() == 4) x[i, , drop=drop] else x[i, drop=drop]
383      else      else
384          x[i,j, drop=drop]          x[i,j, drop=drop]
385      if(isDiagonal(x)) as(x, "diagonalMatrix") else x      if(isS4(x) && isDiagonal(x)) as(x, "diagonalMatrix") else x
386  }  }
387    
388  setMethod("[", signature(x = "diagonalMatrix", i = "index",  setMethod("[", signature(x = "diagonalMatrix", i = "index",
389                           j = "index", drop = "logical"), subDiag)                           j = "index", drop = "logical"), subDiag)
390  setMethod("[", signature(x = "diagonalMatrix", i = "index",  setMethod("[", signature(x = "diagonalMatrix", i = "index",
391                          j = "missing", drop = "logical"),                          j = "missing", drop = "logical"),
392            function(x, i, drop) subDiag(x, i=i, drop=drop))            function(x, i, j, ..., drop) {
393                  na <- nargs()
394                  Matrix.msg("diag[i,m,l] : nargs()=", na, .M.level = 2)
395                  if(na == 4)
396                       subDiag(x, i=i, , drop=drop)
397                  else subDiag(x, i=i,   drop=drop)
398              })
399  setMethod("[", signature(x = "diagonalMatrix", i = "missing",  setMethod("[", signature(x = "diagonalMatrix", i = "missing",
400                           j = "index", drop = "logical"),                           j = "index", drop = "logical"),
401            function(x, j, drop) subDiag(x, j=j, drop=drop))            function(x, i, j, ..., drop) subDiag(x, j=j, drop=drop))
402    
403  ## When you assign to a diagonalMatrix, the result should be  ## When you assign to a diagonalMatrix, the result should be
404  ## diagonal or sparse ---  ## diagonal or sparse ---
405  ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch  ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch
406  replDiag <- function(x, i, j, value) {  ## Only(?) current bug:  x[i] <- value  is wrong when  i is *vector*
407      x <- as(x, "sparseMatrix")  replDiag <- function(x, i, j, ..., value) {
408        x <- as(x, "CsparseMatrix")# was "Tsparse.." till 2012-07
409      if(missing(i))      if(missing(i))
410          x[, j] <- value          x[, j] <- value
411      else if(missing(j))      else if(missing(j)) { ##  x[i , ] <- v  *OR*   x[i] <- v
412            na <- nargs()
413    ##         message("diagnosing replDiag() -- nargs()= ", na)
414            if(na == 4)
415          x[i, ] <- value          x[i, ] <- value
416      else          else if(na == 3)
417                x[i] <- value
418            else stop("Internal bug: nargs()=",na,"; please report")
419        } else
420          x[i,j] <- value          x[i,j] <- value
421      if(isDiagonal(x)) as(x, "diagonalMatrix") else x      if(isDiagonal(x)) as(x, "diagonalMatrix") else x
422  }  }
423    
424  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index",
425                                  j = "index", value = "replValue"), replDiag)                                  j = "index", value = "replValue"), replDiag)
426    
427  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index",
428                                  j = "missing", value = "replValue"),                                  j = "missing", value = "replValue"),
429                   function(x, i, value) replDiag(x, i=i, value=value))                   function(x,i,j, ..., value) {
430                         ## message("before replDiag() -- nargs()= ", nargs())
431                         if(nargs() == 3)
432                             replDiag(x, i=i, value=value)
433                         else ## nargs() == 4 :
434                             replDiag(x, i=i, , value=value)
435                     })
436    
437    setReplaceMethod("[", signature(x = "diagonalMatrix",
438                                    i = "matrix", # 2-col.matrix
439                                    j = "missing", value = "replValue"),
440                     function(x,i,j, ..., value) {
441                         if(ncol(i) == 2) {
442                             if(all((ii <- i[,1]) == i[,2])) { # replace in diagonal only
443                                 if(x@diag == "U") {
444                                     one <- as1(x@x)
445                                     if(any(value != one | is.na(value))) {
446                                         x@diag <- "N"
447                                         x@x <- rep.int(one, x@Dim[1])
448                                     } else return(x)
449                                 }
450                                 x@x[ii] <- value
451                                 x
452                             } else { ## no longer diagonal, but remain sparse:
453                                 x <- as(x, "TsparseMatrix")
454                                 x[i] <- value
455                                 x
456                             }
457                         }
458                         else { # behave as "base R": use as if vector
459                             x <- as(x, "matrix")
460                             x[i] <- value
461                             Matrix(x)
462                         }
463                     })
464    
465  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing",
466                                  j = "index", value = "replValue"),                                  j = "index", value = "replValue"),
467                   function(x, j, value) replDiag(x, j=j, value=value))                   function(x,i,j, ..., value) replDiag(x, j=j, value=value))
468    
469    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",
470                                    value = "sparseMatrix"),
471                     function (x, i, j, ..., value)
472                     callGeneric(x=x, , j=j, value = as(value, "sparseVector")))
473    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "missing",
474                                    value = "sparseMatrix"),
475                     function (x, i, j, ..., value)
476                     callGeneric(x=x, i=i, , value = as(value, "sparseVector")))
477    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "index",
478                                    value = "sparseMatrix"),
479                     function (x, i, j, ..., value)
480                     callGeneric(x=x, i=i, j=j, value = as(value, "sparseVector")))
481    
482    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",
483                                    value = "sparseVector"),
484                     replDiag)
485    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "missing",
486                                    value = "sparseVector"),
487                     replDiag)
488    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "index",
489                                    value = "sparseVector"),
490                     replDiag)
491    
492    
493  setMethod("t", signature(x = "diagonalMatrix"),  setMethod("t", signature(x = "diagonalMatrix"),
# Line 232  Line 498 
498  setMethod("isTriangular", signature(object = "diagonalMatrix"),  setMethod("isTriangular", signature(object = "diagonalMatrix"),
499            function(object) TRUE)            function(object) TRUE)
500  setMethod("isSymmetric", signature(object = "diagonalMatrix"),  setMethod("isSymmetric", signature(object = "diagonalMatrix"),
501            function(object) TRUE)            function(object, ...) TRUE)
502    
503    setMethod("symmpart", signature(x = "diagonalMatrix"), function(x) x)
504    setMethod("skewpart", signature(x = "diagonalMatrix"), setZero)
505    
506  setMethod("chol", signature(x = "ddiMatrix"),# pivot = "ANY"  setMethod("chol", signature(x = "ddiMatrix"),
507            function(x, pivot) {            function(x, pivot, ...) {
508                if(any(x@x < 0)) stop("chol() is undefined for diagonal matrix with negative entries")                if(x@diag == "U") return(x)
509                  ## else
510                  if(any(x@x < 0))
511                      stop("chol() is undefined for diagonal matrix with negative entries")
512                x@x <- sqrt(x@x)                x@x <- sqrt(x@x)
513                x                x
514            })            })
515  ## chol(L) is L for logical diagonal:  ## chol(L) is L for logical diagonal:
516  setMethod("chol", signature(x = "ldiMatrix"), function(x, pivot) x)  setMethod("chol", signature(x = "ldiMatrix"), function(x, pivot, ...) x)
517    
518  setMethod("!", "ldiMatrix", function(e1) {  setMethod("determinant", signature(x = "diagonalMatrix", logarithm = "logical"),
519      if(e1@diag == "N")            function(x, logarithm, ...) mkDet(.diag.x(x), logarithm))
520          e1@x <- !e1@x  
521      else { ## "U"  setMethod("norm", signature(x = "diagonalMatrix", type = "character"),
522          e1@diag <- "N"            function(x, type, ...) {
523          e1@x <- rep.int(FALSE, e1@Dim[1])                if((n <- x@Dim[1]) == 0) return(0) # as for "sparseMatrix"
524                  type <- toupper(substr(type[1], 1, 1))
525                  isU <- (x@diag == "U") # unit-diagonal
526                  if(type == "F") sqrt(if(isU) n else sum(x@x^2))
527                  else { ## norm == "I","1","O","M" :
528                      if(isU) 1 else max(abs(x@x))
529      }      }
     e1  
530  })  })
531    
532    
533    
534  ## Basic Matrix Multiplication {many more to add}  ## Basic Matrix Multiplication {many more to add}
535  ##       ---------------------  ##       ---------------------
536  ## Note that "ldi" logical are treated as numeric  ## Note that "ldi" logical are treated as numeric
537  diagdiagprod <- function(x, y) {  diagdiagprod <- function(x, y) {
538      if(any(dim(x) != dim(y))) stop("non-matching dimensions")      dimCheck(x,y)
539      if(x@diag != "U") {      if(x@diag != "U") {
540          if(y@diag != "U") {          if(y@diag != "U") {
541              nx <- x@x * y@x              nx <- x@x * y@x
# Line 285  Line 563 
563    
564    
565  diagmatprod <- function(x, y) {  diagmatprod <- function(x, y) {
566        ## x is diagonalMatrix
567      dx <- dim(x)      dx <- dim(x)
568      dy <- dim(y)      dy <- dim(y)
569      if(dx[2] != dy[1]) stop("non-matching dimensions")      if(dx[2] != dy[1]) stop("non-matching dimensions")
     n <- dx[1]  
570      as(if(x@diag == "U") y else x@x * y, "Matrix")      as(if(x@diag == "U") y else x@x * y, "Matrix")
571  }  }
   
572  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),
573            diagmatprod)            diagmatprod)
574    ## sneaky .. :
575  formals(diagmatprod) <- alist(x=, y=NULL)  formals(diagmatprod) <- alist(x=, y=NULL)
576  setMethod("crossprod", signature(x = "diagonalMatrix", y = "matrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "matrix"),
577            diagmatprod)            diagmatprod)
 setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "matrix"),  
           diagmatprod)  
578    
579  diagdgeprod <- function(x, y) {  diagGeprod <- function(x, y) {
580      dx <- dim(x)      dx <- dim(x)
581      dy <- dim(y)      dy <- dim(y)
582      if(dx[2] != dy[1]) stop("non-matching dimensions")      if(dx[2] != dy[1]) stop("non-matching dimensions")
# Line 308  Line 584 
584          y@x <- x@x * y@x          y@x <- x@x * y@x
585      y      y
586  }  }
587  setMethod("%*%", signature(x = "diagonalMatrix", y = "dgeMatrix"),  setMethod("%*%", signature(x= "diagonalMatrix", y= "dgeMatrix"), diagGeprod)
588            diagdgeprod, valueClass = "dgeMatrix")  setMethod("%*%", signature(x= "diagonalMatrix", y= "lgeMatrix"), diagGeprod)
589  formals(diagdgeprod) <- alist(x=, y=NULL)  formals(diagGeprod) <- alist(x=, y=NULL)
590  setMethod("crossprod", signature(x = "diagonalMatrix", y = "dgeMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "dgeMatrix"),
591            diagdgeprod, valueClass = "dgeMatrix")            diagGeprod, valueClass = "dgeMatrix")
592    setMethod("crossprod", signature(x = "diagonalMatrix", y = "lgeMatrix"),
593              diagGeprod)
594    
595  setMethod("%*%", signature(x = "matrix", y = "diagonalMatrix"),  matdiagprod <- function(x, y) {
           function(x, y) {  
596                dx <- dim(x)                dx <- dim(x)
597                dy <- dim(y)                dy <- dim(y)
598                if(dx[2] != dy[1]) stop("non-matching dimensions")                if(dx[2] != dy[1]) stop("non-matching dimensions")
599                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]))
600            })  }
601    setMethod("%*%", signature(x = "matrix", y = "diagonalMatrix"),
602              matdiagprod)
603    formals(matdiagprod) <- alist(x=, y=NULL)
604    setMethod("tcrossprod", signature(x = "matrix", y = "diagonalMatrix"),
605              matdiagprod)
606    
607  setMethod("%*%", signature(x = "dgeMatrix", y = "diagonalMatrix"),  gediagprod <- function(x, y) {
           function(x, y) {  
608                dx <- dim(x)                dx <- dim(x)
609                dy <- dim(y)                dy <- dim(y)
610                if(dx[2] != dy[1]) stop("non-matching dimensions")                if(dx[2] != dy[1]) stop("non-matching dimensions")
611                if(y@diag == "N")                if(y@diag == "N")
612                    x@x <- x@x * rep(y@x, each = dx[1])                    x@x <- x@x * rep(y@x, each = dx[1])
613                x                x
614            })  }
615    setMethod("%*%", signature(x= "dgeMatrix", y= "diagonalMatrix"), gediagprod)
616    setMethod("%*%", signature(x= "lgeMatrix", y= "diagonalMatrix"), gediagprod)
617    formals(gediagprod) <- alist(x=, y=NULL)
618    setMethod("tcrossprod", signature(x = "dgeMatrix", y = "diagonalMatrix"),
619              gediagprod)
620    setMethod("tcrossprod", signature(x = "lgeMatrix", y = "diagonalMatrix"),
621              gediagprod)
622    
623  ## crossprod {more of these}  ## crossprod {more of these}
624    
625  ## tcrossprod --- all are not yet there: do the dense ones here:  ## tcrossprod --- all are not yet there: do the dense ones here:
626    
627    setMethod("%*%", signature(x = "diagonalMatrix", y = "denseMatrix"),
628              function(x, y) if(x@diag == "U") y else x %*% as(y, "generalMatrix"))
629    setMethod("%*%", signature(x = "denseMatrix", y = "diagonalMatrix"),
630              function(x, y) if(y@diag == "U") x else as(x, "generalMatrix") %*% y)
631    
632    
633  ## FIXME:  ## FIXME:
634  ## setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "denseMatrix"),  ## setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "denseMatrix"),
635  ##        function(x, y = NULL) {  ##        function(x, y = NULL) {
636  ##           })  ##           })
637    
638  ## setMethod("tcrossprod", signature(x = "denseMatrix", y = "diagonalMatrix"),  Cspdiagprod <- function(x, y) {
639  ##        function(x, y = NULL) {      dx <- dim(x <- .Call(Csparse_diagU2N, x))
640  ##           })      dy <- dim(y)
641        if(dx[2] != dy[1]) stop("non-matching dimensions")
642        if(y@diag == "N") {
643            if(!all(y@x[1L] == y@x[-1L]) && is(x, "symmetricMatrix"))
644                x <- as(x, "generalMatrix")
645            ind <- rep.int(seq_len(dx[2]), x@p[-1] - x@p[-dx[2]-1L])
646            x@x <- x@x * y@x[ind]
647        }
648        if(is(x, "compMatrix") && length(xf <- x@factors)) {
649            ## instead of dropping all factors, be smart about some
650            ## TODO ......
651            x@factors <- list()
652        }
653        x
654    }
655    
656    diagCspprod <- function(x, y) {
657        dx <- dim(x)
658        dy <- dim(y <- .Call(Csparse_diagU2N, y))
659        if(dx[2] != dy[1]) stop("non-matching dimensions")
660        if(x@diag == "N") {
661            if(!all(x@x[1L] == x@x[-1L]) && is(y, "symmetricMatrix"))
662                y <- as(y, "generalMatrix")
663            y@x <- y@x * x@x[y@i + 1L]
664        }
665        if(is(y, "compMatrix") && length(yf <- y@factors)) {
666            ## instead of dropping all factors, be smart about some
667            ## TODO
668            keep <- character()
669            if(iLU <- names(yf) == "LU") {
670                ## TODO keep <- "LU"
671            }
672            y@factors <- yf[keep]
673        }
674        y
675    }
676    
677    setMethod("crossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
678              function(x, y = NULL) diagCspprod(x, y))
679    
680  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
681            function(x, y = NULL) { x <- as(x, "sparseMatrix"); callGeneric() })            function(x, y = NULL) diagCspprod(x, as(y, "CsparseMatrix")))
682    
683    ## Prefer calling diagCspprod to Cspdiagprod if going to transpose anyway
684    ##  x'y == (y'x)'
685    setMethod("crossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
686              function(x, y = NULL) t(diagCspprod(y, x)))
687    
688  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
689            function(x, y = NULL) { y <- as(y, "sparseMatrix"); callGeneric() })            function(x, y = NULL) t(diagCspprod(y, as(x, "Csparsematrix"))))
690    
691    setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
692              function(x, y = NULL) diagCspprod(x, t(y)))
693    
694  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
695            function(x, y = NULL) { x <- as(x, "sparseMatrix"); callGeneric() })            function(x, y = NULL) diagCspprod(x, t(as(y, "CsparseMatrix"))))
696    
697    setMethod("tcrossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
698              function(x, y = NULL) Cspdiagprod(x, y))
699    
700  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
701            function(x, y = NULL) { y <- as(y, "sparseMatrix"); callGeneric() })            function(x, y = NULL) Cspdiagprod(as(x, "CsparseMatrix"), y))
702    
703    setMethod("%*%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
704              function(x, y) diagCspprod(x, y))
705    
 ## FIXME?: In theory, this can be done *FASTER*, in some cases, via tapply1()  
706  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
707            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)!  
708    
709  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
710            function(x, y) x %*% as(y, "sparseMatrix"))            function(x, y) Cspdiagprod(as(x, "CsparseMatrix"), y))
711    
712    setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
713              function(x, y) Cspdiagprod(x, y))
714    
715    ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*
716    ##       do indeed work by going through sparse (and *not* ddense)!
717    
718  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),
719            function(a, b, ...) {            function(a, b, ...) {
# Line 384  Line 723 
723            })            })
724    
725  solveDiag <- function(a, b, ...) {  solveDiag <- function(a, b, ...) {
726      if((n <- a@Dim[1]) != nrow(b))      if(a@Dim[1] != nrow(b))
727          stop("incompatible matrix dimensions")          stop("incompatible matrix dimensions")
728      ## trivially invert a 'in place' and multiply:      ## trivially invert a 'in place' and multiply:
729      a@x <- 1/ a@x      a@x <- 1/ a@x
# Line 396  Line 735 
735  setMethod("solve", signature(a = "diagonalMatrix", b = "Matrix"),  setMethod("solve", signature(a = "diagonalMatrix", b = "Matrix"),
736            solveDiag)            solveDiag)
737    
738    ## Schur()  ---> ./eigen.R
739    
740    
741    
742  ### ---------------- diagonal  o  sparse  -----------------------------  ###---------------- <Ops> (<Arith>, <Logic>, <Compare> ) ----------------------
   
743    
744  ## Use function for several signatures, in order to evade  ## Use function for several signatures, in order to evade
745  ## ambiguous dispatch for "ddi", since there's also Arith(ddense., ddense.)  diagOdiag <- function(e1,e2) {
746  diagOdiag <- function(e1,e2) { # result should also be diagonal      ## result should also be diagonal _ if possible _
747      r <- callGeneric(.diag.x(e1), .diag.x(e2)) # error if not "compatible"      r <- callGeneric(.diag.x(e1), .diag.x(e2)) # error if not "compatible"
748      if(is.numeric(r)) {      ## Check what happens with non-diagonals, i.e. (0 o 0), (FALSE o 0), ...:
749        r00 <- callGeneric(if(is.numeric(e1@x)) 0 else FALSE,
750                           if(is.numeric(e2@x)) 0 else FALSE)
751        if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* diagonal
752            if(is.numeric(r)) { # "double" *or* "integer"
753          if(is.numeric(e2@x)) {          if(is.numeric(e2@x)) {
754              e2@x <- r; return(.diag.2N(e2)) }              e2@x <- r; return(.diag.2N(e2)) }
755          if(!is.numeric(e1@x))          if(!is.numeric(e1@x))
756              ## e.g. e1, e2 are logical;              ## e.g. e1, e2 are logical;
757              e1 <- as(e1, "dMatrix")              e1 <- as(e1, "dMatrix")
758                if(!is.double(r)) r <- as.double(r)
759      }      }
760      else if(is.logical(r))      else if(is.logical(r))
761          e1 <- as(e1, "lMatrix")          e1 <- as(e1, "lMatrix")
# Line 419  Line 763 
763      e1@x <- r      e1@x <- r
764      .diag.2N(e1)      .diag.2N(e1)
765  }  }
766        else { ## result not diagonal, but at least symmetric:
767            ## e.g., m == m
768            isNum <- (is.numeric(r) || is.numeric(r00))
769            isLog <- (is.logical(r) || is.logical(r00))
770            Matrix.msg("exploding <diag> o <diag> into dense matrix", .M.level = 2)
771            d <- e1@Dim
772            n <- d[1]
773            stopifnot(length(r) == n)
774            if(isNum && !is.double(r)) r <- as.double(r)
775            xx <- as.vector(matrix(rbind(r, matrix(r00,n,n)), n,n))
776            newcl <-
777                paste0(if(isNum) "d" else if(isLog) {
778                    if(!any(is.na(r)) && !any(is.na(r00))) "n" else "l"
779                } else stop("not yet implemented .. please report"), "syMatrix")
780    
781            new(newcl, Dim = e1@Dim, Dimnames = e1@Dimnames, x = xx)
782        }
783    }
784    
785    ### This would be *the* way, but we get tons of "ambiguous method dispatch"
786    ## we use this hack instead of signature  x = "diagonalMatrix" :
787    diCls <- names(getClass("diagonalMatrix")@subclasses)
788    if(FALSE) {
789  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "diagonalMatrix"),  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "diagonalMatrix"),
790            diagOdiag)            diagOdiag)
791  ## These two are just for method disambiguation:  } else { ## These are just for method disambiguation:
792  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "diagonalMatrix"),      for(c1 in diCls)
793            diagOdiag)          for(c2 in diCls)
794  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "ddiMatrix"),              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)
795            diagOdiag)  }
796    
797    ## diagonal  o  triangular  |-->  triangular
798    ## diagonal  o  symmetric   |-->  symmetric
799    ##    {also when other is sparse: do these "here" --
800    ##     before conversion to sparse, since that loses "diagonality"}
801    diagOtri <- function(e1,e2) {
802        ## result must be triangular
803        r <- callGeneric(d1 <- .diag.x(e1), diag(e2)) # error if not "compatible"
804        ## Check what happens with non-diagonals, i.e. (0 o 0), (FALSE o 0), ...:
805        e1.0 <- if(.n1 <- is.numeric(d1   )) 0 else FALSE
806        r00 <- callGeneric(e1.0, if(.n2 <- is.numeric(e2[0])) 0 else FALSE)
807        if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* triangular
808            diag(e2) <- r
809            ## check what happens "in the triangle"
810            e2.2 <- if(.n2) 2 else TRUE
811            if(!callGeneric(e1.0, e2.2) == e2.2) { # values "in triangle" can change:
812                n <- dim(e2)[1L]
813                it <- indTri(n, upper = (e2@uplo == "U"))
814                e2[it] <- callGeneric(e1.0, e2[it])
815            }
816            e2
817        }
818        else { ## result not triangular ---> general
819            rr <- as(e2, "generalMatrix")
820            diag(rr) <- r
821            rr
822        }
823    }
824    
825    
826    setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "triangularMatrix"),
827              diagOtri)
828    ## For the reverse,  Ops == "Arith" | "Compare" | "Logic"
829    ##   'Arith'  :=  '"+"', '"-"', '"*"', '"^"', '"%%"', '"%/%"', '"/"'
830    setMethod("Arith", signature(e1 = "triangularMatrix", e2 = "diagonalMatrix"),
831              function(e1,e2)
832          { ## this must only trigger for *dense* e1
833              switch(.Generic,
834                     "+" = .Call(dtrMatrix_addDiag, as(e1,"dtrMatrix"),   .diag.x(e2)),
835                     "-" = .Call(dtrMatrix_addDiag, as(e1,"dtrMatrix"), - .diag.x(e2)),
836                     "*" = {
837                         n <- e2@Dim[1L]
838                         d2 <- if(e2@diag == "U") { # unit-diagonal
839                             d <- rep.int(as1(e2@x), n)
840                             e2@x <- d
841                             e2@diag <- "N"
842                             d
843                         } else e2@x
844                         e2@x <- diag(e1) * d2
845                         e2
846                     },
847                     "^" = { ## will be dense ( as  <ANY> ^ 0 == 1 ):
848                         e1 ^ as(e2, "denseMatrix")
849                     },
850                     ## otherwise:
851                     callGeneric(e1, diag2Tsmart(e2,e1)))
852    })
853    
854    ## Compare --> 'swap' (e.g.   e1 < e2   <==>  e2 > e1 ):
855    setMethod("Compare", signature(e1 = "triangularMatrix", e2 = "diagonalMatrix"),
856              .Cmp.swap)
857    ## '&' and "|'  are commutative:
858    setMethod("Logic", signature(e1 = "triangularMatrix", e2 = "diagonalMatrix"),
859              function(e1,e2) callGeneric(e2, e1))
860    
861  ## For almost everything else, diag* shall be treated "as sparse" :  ## For almost everything else, diag* shall be treated "as sparse" :
862  ## These are cheap implementations via coercion  ## These are cheap implementations via coercion
863    
864  ## for disambiguation  ## For disambiguation --- define this for "sparseMatrix" , then for "ANY";
865  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "sparseMatrix"),  ## and because we can save an .M.kind() call, we use this explicit
866            function(e1,e2) callGeneric(as(e1, "sparseMatrix"), e2))  ## "hack" for all diagonalMatrix *subclasses* instead of just "diagonalMatrix" :
867  setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "diagonalMatrix"),  ##
868            function(e1,e2) callGeneric(e1, as(e2, "sparseMatrix")))  ## ddi*:
869  ## in general:  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "sparseMatrix"),
870  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "ANY"),            function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2))
871            function(e1,e2) callGeneric(as(e1,"sparseMatrix"), e2))  setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "ddiMatrix"),
872  setMethod("Ops", signature(e1 = "ANY", e2 = "diagonalMatrix"),            function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "d")))
873            function(e1,e2) callGeneric(e1, as(e2,"sparseMatrix")))  ## ldi*
874    setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "sparseMatrix"),
875              function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2))
876    setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "ldiMatrix"),
877              function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "l")))
878    
879    ## Ops:  Arith  --> numeric : "dMatrix"
880    ##       Compare --> logical
881    ##       Logic   --> logical: "lMatrix"
882    
883    ## Other = "numeric" : stay diagonal if possible
884    ## ddi*: Arith: result numeric, potentially ddiMatrix
885    for(arg2 in c("numeric","logical"))
886    setMethod("Arith", signature(e1 = "ddiMatrix", e2 = arg2),
887              function(e1,e2) {
888                  n <- e1@Dim[1]
889                  f0 <- callGeneric(0, e2)
890                  if(all(is0(f0))) { # remain diagonal
891                      L1 <- (le <- length(e2)) == 1L
892                      if(e1@diag == "U") {
893                          if(any((r <- callGeneric(1, e2)) != 1)) {
894                              e1@diag <- "N"
895                              e1@x[seq_len(n)] <- r # possibly recycling r
896                          } ## else: result = e1  (is "U" diag)
897                      } else {
898                          r <- callGeneric(e1@x, e2)
899                          ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
900                          e1@x[] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
901                      }
902                      e1
903                  } else
904                      callGeneric(diag2tT.u(e1,e2, "d"), e2)
905              })
906    
907    for(arg1 in c("numeric","logical"))
908    setMethod("Arith", signature(e1 = arg1, e2 = "ddiMatrix"),
909              function(e1,e2) {
910                  n <- e2@Dim[1]
911                  f0 <- callGeneric(e1, 0)
912                  if(all(is0(f0))) { # remain diagonal
913                      L1 <- (le <- length(e1)) == 1L
914                      if(e2@diag == "U") {
915                          if(any((r <- callGeneric(e1, 1)) != 1)) {
916                              e2@diag <- "N"
917                              e2@x[seq_len(n)] <- r # possibly recycling r
918                          } ## else: result = e2  (is "U" diag)
919                      } else {
920                          r <- callGeneric(e1, e2@x)
921                          ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
922                          e2@x[] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
923                      }
924                      e2
925                  } else
926                      callGeneric(e1, diag2tT.u(e2,e1, "d"))
927              })
928    
929    ## ldi* Arith --> result numeric, potentially ddiMatrix
930    for(arg2 in c("numeric","logical"))
931    setMethod("Arith", signature(e1 = "ldiMatrix", e2 = arg2),
932              function(e1,e2) {
933                  n <- e1@Dim[1]
934                  f0 <- callGeneric(0, e2)
935                  if(all(is0(f0))) { # remain diagonal
936                      L1 <- (le <- length(e2)) == 1L
937                      E <- copyClass(e1, "ddiMatrix", c("diag", "Dim", "Dimnames"))#FIXME: if ok, check=FALSE
938                      ## E <- copyClass(e1, "ddiMatrix", check=FALSE)
939                      ## storage.mode(E@x) <- "double"
940                      if(e1@diag == "U") {
941                          if(any((r <- callGeneric(1, e2)) != 1)) {
942                              E@diag <- "N"
943                              E@x[seq_len(n)] <- r # possibly recycling r
944                          } ## else: result = E  (is "U" diag)
945                      } else {
946                          r <- callGeneric(e1@x, e2)
947                          ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
948                          E@x[seq_len(n)] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
949                      }
950                      E
951                  } else
952                      callGeneric(diag2tT.u(e1,e2, "l"), e2)
953              })
954    
955    for(arg1 in c("numeric","logical"))
956    setMethod("Arith", signature(e1 = arg1, e2 = "ldiMatrix"),
957              function(e1,e2) {
958                  n <- e2@Dim[1]
959                  f0 <- callGeneric(e1, 0)
960                  if(all(is0(f0))) { # remain diagonal
961                      L1 <- (le <- length(e1)) == 1L
962                      E <- copyClass(e2, "ddiMatrix", c("diag", "Dim", "Dimnames"))#FIXME: if ok, check=FALSE
963                      ## E <- copyClass(e2, "ddiMatrix", check=FALSE)
964                      ## storage.mode(E@x) <- "double"
965                      if(e2@diag == "U") {
966                          if(any((r <- callGeneric(e1, 1)) != 1)) {
967                              E@diag <- "N"
968                              E@x[seq_len(n)] <- r # possibly recycling r
969                          } ## else: result = E  (is "U" diag)
970                      } else {
971                          r <- callGeneric(e1, e2@x)
972                          ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
973                          E@x[seq_len(n)] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
974                      }
975                      E
976                  } else
977                      callGeneric(e1, diag2tT.u(e2,e1, "l"))
978              })
979    
980    ## ddi*: for "Ops" without "Arith": <Compare> or <Logic> --> result logical, potentially ldi
981    ##
982    ## Note that  ("numeric", "ddiMatrix")  is simply swapped, e.g.,
983    if(FALSE) {
984        selectMethod("<", c("numeric","lMatrix"))# Compare
985        selectMethod("&", c("numeric","lMatrix"))# Logic
986    } ## so we don't need to define a method here :
987    
988    for(arg2 in c("numeric","logical"))
989    setMethod("Ops", signature(e1 = "ddiMatrix", e2 = arg2),
990              function(e1,e2) {
991                  n <- e1@Dim[1]
992                  f0 <- callGeneric(0, e2)
993                  if(all(is0(f0))) { # remain diagonal
994                      L1 <- (le <- length(e2)) == 1L
995                      E <- copyClass(e1, "ldiMatrix", c("diag", "Dim", "Dimnames"))#FIXME: if ok, check=FALSE
996                      ## E <- copyClass(e1, "ldiMatrix", check=FALSE)
997                      ## storage.mode(E@x) <- "logical"
998                      if(e1@diag == "U") {
999                          if(any((r <- callGeneric(1, e2)) != 1)) {
1000                              E@diag <- "N"
1001                              E@x[seq_len(n)] <- r # possibly recycling r
1002                          } ## else: result = E  (is "U" diag)
1003                      } else {
1004                          r <- callGeneric(e1@x, e2)
1005                          ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
1006                          E@x[seq_len(n)] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
1007                      }
1008                      E
1009                  } else
1010                      callGeneric(diag2tT.u(e1,e2, "d"), e2)
1011              })
1012    
1013    ## ldi*: for "Ops" without "Arith": <Compare> or <Logic> --> result logical, potentially ldi
1014    for(arg2 in c("numeric","logical"))
1015    setMethod("Ops", signature(e1 = "ldiMatrix", e2 = arg2),
1016              function(e1,e2) {
1017                  n <- e1@Dim[1]
1018                  f0 <- callGeneric(FALSE, e2)
1019                  if(all(is0(f0))) { # remain diagonal
1020                      L1 <- (le <- length(e2)) == 1L
1021    
1022                      if(e1@diag == "U") {
1023                          if(any((r <- callGeneric(TRUE, e2)) != 1)) {
1024                              e1@diag <- "N"
1025                              e1@x[seq_len(n)] <- r # possibly recycling r
1026                          } ## else: result = e1  (is "U" diag)
1027                      } else {
1028                          r <- callGeneric(e1@x, e2)
1029                          ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
1030                          e1@x[] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
1031                      }
1032                      e1
1033                  } else
1034                      callGeneric(diag2tT.u(e1,e2, "l"), e2)
1035              })
1036    
1037    
1038    ## Not {"sparseMatrix", "numeric} :  {"denseMatrix", "matrix", ... }
1039    for(other in c("ANY", "Matrix", "dMatrix")) {
1040        ## ddi*:
1041        setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other),
1042                  function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2))
1043        setMethod("Ops", signature(e1 = other, e2 = "ddiMatrix"),
1044                  function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "d")))
1045        ## ldi*:
1046        setMethod("Ops", signature(e1 = "ldiMatrix", e2 = other),
1047                  function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2))
1048        setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),
1049                  function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "l")))
1050    }
1051    
1052    ## Direct subclasses of "denseMatrix": currently ddenseMatrix, ldense... :
1053    dense.subCl <- local({ dM.scl <- getClass("denseMatrix")@subclasses
1054                           names(dM.scl)[vapply(dM.scl, slot, 0, "distance") == 1] })
1055    for(DI in diCls) {
1056        dMeth <- if(extends(DI, "dMatrix"))
1057            function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2)
1058        else # "lMatrix", the only other kind for now
1059            function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2)
1060        for(c2 in c(dense.subCl, "Matrix")) {
1061            for(Fun in c("*", "&")) {
1062                setMethod(Fun, signature(e1 = DI, e2 = c2),
1063                          function(e1,e2) callGeneric(e1, Diagonal(x = diag(e2))))
1064                setMethod(Fun, signature(e1 = c2, e2 = DI),
1065                          function(e1,e2) callGeneric(Diagonal(x = diag(e1)), e2))
1066            }
1067            setMethod("^", signature(e1 = c2, e2 = DI),
1068                      function(e1,e2) callGeneric(Diagonal(x = diag(e1)), e2))
1069            for(Fun in c("%%", "%/%", "/")) ## 0 <op> 0 |--> NaN  for these.
1070                setMethod(Fun, signature(e1 = DI, e2 = c2), dMeth)
1071        }
1072    }
1073    
1074    
1075    ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"
1076    ### ----------   the last 4: separately here
1077    for(cl in diCls) {
1078    setMethod("any", cl,
1079              function (x, ..., na.rm) {
1080                  if(any(x@Dim == 0)) FALSE
1081                  else if(x@diag == "U") TRUE else any(x@x, ..., na.rm = na.rm)
1082              })
1083    setMethod("all",  cl, function (x, ..., na.rm) {
1084        n <- x@Dim[1]
1085        if(n >= 2) FALSE
1086        else if(n == 0 || x@diag == "U") TRUE
1087        else all(x@x, ..., na.rm = na.rm)
1088    })
1089    setMethod("prod", cl, function (x, ..., na.rm) {
1090        n <- x@Dim[1]
1091        if(n >= 2) 0
1092        else if(n == 0 || x@diag == "U") 1
1093        else ## n == 1, diag = "N" :
1094            prod(x@x, ..., na.rm = na.rm)
1095    })
1096    
1097    setMethod("sum", cl,
1098              function(x, ..., na.rm) {
1099                  r <- sum(x@x, ..., na.rm = na.rm)# double or integer, correctly
1100                  if(x@diag == "U" && !is.na(r)) r + x@Dim[1] else r
1101              })
1102    }
1103    
1104    ## The remaining ones are  max, min, range :
1105    
1106    setMethod("Summary", "ddiMatrix",
1107              function(x, ..., na.rm) {
1108                  if(any(x@Dim == 0)) callGeneric(numeric(0), ..., na.rm=na.rm)
1109                  else if(x@diag == "U")
1110                      callGeneric(x@x, 0, 1, ..., na.rm=na.rm)
1111                  else callGeneric(x@x, 0, ..., na.rm=na.rm)
1112              })
1113    setMethod("Summary", "ldiMatrix",
1114              function(x, ..., na.rm) {
1115                  if(any(x@Dim == 0)) callGeneric(logical(0), ..., na.rm=na.rm)
1116                  else if(x@diag == "U")
1117                      callGeneric(x@x, FALSE, TRUE, ..., na.rm=na.rm)
1118                  else callGeneric(x@x, FALSE, ..., na.rm=na.rm)
1119              })
1120    
1121    
1122    
# Line 455  Line 1131 
1131      invisible(x)      invisible(x)
1132  }  }
1133    
1134    ## somewhat consistent with "print" for sparseMatrix :
1135    setMethod("print", signature(x = "diagonalMatrix"), prDiag)
1136    
1137  setMethod("show", signature(object = "diagonalMatrix"),  setMethod("show", signature(object = "diagonalMatrix"),
1138            function(object) {            function(object) {
1139                d <- dim(object)                d <- dim(object)
1140                cl <- class(object)                cl <- class(object)
1141                cat(sprintf('%d x %d diagonal matrix of class "%s"\n',                cat(sprintf('%d x %d diagonal matrix of class "%s"',
1142                            d[1], d[2], cl))                            d[1], d[2], cl))
1143                  if(d[1] < 50) {
1144                      cat("\n")
1145                prDiag(object)                prDiag(object)
1146                  } else {
1147                      cat(", with diagonal entries\n")
1148                      show(diag(object))
1149                      invisible(object)
1150                  }
1151            })            })
1152    
1153    rm(dense.subCl, diCls)# not used elsewhere
1154    
1155    setMethod("summary", signature(object = "diagonalMatrix"),
1156              function(object, ...) {
1157                  d <- dim(object)
1158                  r <- summary(object@x, ...)
1159                  attr(r, "header") <-
1160                      sprintf('%d x %d diagonal Matrix of class "%s"',
1161                              d[1], d[2], class(object))
1162                  ## use ole' S3 technology for such a simple case
1163                  class(r) <- c("diagSummary", class(r))
1164                  r
1165              })
1166    
1167    print.diagSummary <- function (x, ...) {
1168        cat(attr(x, "header"),"\n")
1169        class(x) <- class(x)[-1]
1170        print(x, ...)
1171        invisible(x)
1172    }

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

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