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

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

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