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

Legend:
Removed from v.2096  
changed lines
  Added in v.3015

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