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 2764, Fri Feb 17 19:58:29 2012 UTC
# Line 16  Line 16 
16      if(missing(x)) ## unit diagonal matrix      if(missing(x)) ## unit diagonal matrix
17          new("ddiMatrix", Dim = c(n,n), diag = "U")          new("ddiMatrix", Dim = c(n,n), diag = "U")
18      else {      else {
19          stopifnot(length(x) == n)          lx <- length(x)
20            stopifnot(lx == 1 || lx == n) # but keep 'x' short for now
21          if(is.logical(x))          if(is.logical(x))
22              cl <- "ldiMatrix"              cl <- "ldiMatrix"
23          else if(is.numeric(x)) {          else if(is.numeric(x)) {
# Line 26  Line 27 
27          else if(is.complex(x)) {          else if(is.complex(x)) {
28              cl <- "zdiMatrix"  # will not yet work              cl <- "zdiMatrix"  # will not yet work
29          } else stop("'x' has invalid data type")          } else stop("'x' has invalid data type")
30          new(cl, Dim = c(n,n), diag = "N", x = x)          new(cl, Dim = c(n,n), diag = "N",
31                x = if(lx == 1) rep.int(x,n) else x)
32      }      }
33  }  }
34    
35  ### This is modified from a post of Bert Gunter to R-help on  1 Sep 2005.  .sparseDiagonal <- function(n, x = rep.int(1,m), uplo = "U",
36  ### Bert's code built on a post by Andy Liaw who most probably was influenced                              shape = if(missing(cols)) "t" else "g",
37  ### by earlier posts, notably one by Scott Chasalow on S-news, 16 Jan 2002                              kind, cols = if(n) 0:(n - 1L) else integer(0))
38  ### who posted his bdiag() function written in December 1995.  {
39        stopifnot(n == (n. <- as.integer(n)), (n <- n.) >= 0)
40  bdiag <- function(...) {      if(!missing(cols))
41      if(nargs() == 0) return(new("dgCMatrix"))          stopifnot(0 <= (cols <- as.integer(cols)), cols < n)
42      ## else :      m <- length(cols)
43      mlist <- if (nargs() == 1) as.list(...) else list(...)      if(missing(kind))
44      dims <- sapply(mlist, dim)          kind <-
45                if(is.double(x)) "d"
46                else if(is.logical(x)) "l"
47                else { ## for now
48                    storage.mode(x) <- "double"
49                    "d"
50                }
51        else stopifnot(any(kind == c("d","l","n")))
52        if(kind != "n") {
53            if((lx <- length(x)) == 1) x <- rep.int(x, m)
54            else if(lx != m) stop("length(x) must be either 1 or #{cols}")
55        }
56        stopifnot(is.character(shape), nchar(shape) == 1,
57                  any(shape == c("t","s","g"))) # triangular / symmetric / general
58        if(kind == "n") {
59            if(shape == "g")
60                new("ngCMatrix", Dim = c(n,m), i = cols, p = 0:m)
61            else new(paste0("n", shape, "CMatrix"), Dim = c(n,m), uplo = uplo,
62                     i = cols, p = 0:m)
63        }
64        ## kind != "n" -- have x slot :
65        else if(shape == "g")
66            new(paste0(kind, "gCMatrix"), Dim = c(n,m),
67                x = x, i = cols, p = 0:m)
68        else new(paste0(kind, shape, "CMatrix"), Dim = c(n,m), uplo = uplo,
69                 x = x, i = cols, p = 0:m)
70    }
71    
72    ## Pkg 'spdep' had (relatively slow) versions of this as_dsCMatrix_I()
73    .symDiagonal <- function(n, x = rep.int(1,n), uplo = "U")
74        .sparseDiagonal(n, x, uplo, shape = "s")
75    
76    # instead of   diagU2N(as(Diagonal(n), "CsparseMatrix")), diag = "N" in any case:
77    .trDiagonal <- function(n, x = rep.int(1,n), uplo = "U")
78        .sparseDiagonal(n, x, uplo, shape = "t")
79    
80    
81    ## This is modified from a post of Bert Gunter to R-help on  1 Sep 2005.
82    ## Bert's code built on a post by Andy Liaw who most probably was influenced
83    ## by earlier posts, notably one by Scott Chasalow on S-news, 16 Jan 2002
84    ## who posted his bdiag() function written in December 1995.
85    if(FALSE)##--- no longer used:
86    .bdiag <- function(lst) {
87        ## block-diagonal matrix [a dgTMatrix] from list of matrices
88        stopifnot(is.list(lst), length(lst) >= 1)
89        dims <- sapply(lst, dim, USE.NAMES=FALSE)
90      ## make sure we had all matrices:      ## make sure we had all matrices:
91      if(!(is.matrix(dims) && nrow(dims) == 2))      if(!(is.matrix(dims) && nrow(dims) == 2))
92          stop("some arguments are not matrices")          stop("some arguments are not matrices")
93      csdim <- rbind(rep.int(0:0, 2),      csdim <- rbind(rep.int(0L, 2),
94                     apply(sapply(mlist, dim), 1, cumsum))                     apply(dims, 1, cumsum))
95      ret <- new("dgTMatrix", Dim = as.integer(csdim[nrow(csdim),]))      r <- new("dgTMatrix")
96        r@Dim <- as.integer(csdim[nrow(csdim),])
97      add1 <- matrix(1:0, 2,2)      add1 <- matrix(1:0, 2,2)
98      for(i in seq_along(mlist)) {      for(i in seq_along(lst)) {
99          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])
100          if(is.null(dim(indx))) ## non-square matrix          if(is.null(dim(indx))) ## non-square matrix
101              ret[indx[[1]],indx[[2]]] <- mlist[[i]]              r[indx[[1]],indx[[2]]] <- lst[[i]]
102          else ## square matrix          else ## square matrix
103              ret[indx[,1],indx[,2]] <- mlist[[i]]              r[indx[,1], indx[,2]] <- lst[[i]]
104        }
105        r
106    }
107    ## expand(<mer>) needed something like bdiag() for lower-triangular
108    ## (Tsparse) Matrices; hence Doug Bates provided a much more efficient
109    ##  implementation for those; now extended and generalized:
110    .bdiag <- function(lst) {
111        ## block-diagonal matrix [a dgTMatrix] from list of matrices
112        stopifnot(is.list(lst), (nl <- length(lst)) >= 1)
113    
114        Tlst <- lapply(lapply(lst, Matrix:::as_Csp2), # includes "diagU2N"
115                       as, "TsparseMatrix")
116        if(nl == 1) return(Tlst[[1]])
117        ## else
118        i_off <- c(0L, cumsum(sapply(Tlst, nrow)))
119        j_off <- c(0L, cumsum(sapply(Tlst, ncol)))
120    
121        clss <- sapply(Tlst, class)
122        typ <- substr(clss, 2, 2)
123        knd <- substr(clss, 1, 1)
124        sym <- typ == "s" # symmetric ones
125        tri <- typ == "t" # triangular ones
126        use.n <- any(is.n <- knd == "n")
127        if(use.n && !(use.n <- all(is.n))) {
128            Tlst[is.n] <- lapply(Tlst[is.n], as, "lMatrix")
129            knd [is.n] <- "l"
130        }
131        use.l <- !use.n && all(knd == "l")
132        if(all(sym)) { ## result should be *symmetric*
133            uplos <- sapply(Tlst, slot, "uplo") ## either "U" or "L"
134            tLU <- table(uplos)# of length 1 or 2 ..
135            if(length(tLU) == 1) { ## all "U" or all "L"
136                useU <- uplos[1] == "U"
137            } else { ## length(tLU) == 2, counting "L" and "U"
138                useU <- diff(tLU) >= 0
139                if(useU && (hasL <- tLU[1] > 0))
140                    Tlst[hasL] <- lapply(Tlst[hasL], t)
141                else if(!useU && (hasU <- tLU[2] > 0))
142                    Tlst[hasU] <- lapply(Tlst[hasU], t)
143            }
144            if(use.n) { ## return nsparseMatrix :
145                r <- new("nsTMatrix")
146            } else {
147                r <- new(paste0(if(use.l) "l" else "d", "sTMatrix"))
148                r@x <- unlist(lapply(Tlst, slot, "x"))
149            }
150            r@uplo <- if(useU) "U" else "L"
151      }      }
152      ## slightly debatable if we really should return Csparse.. :      else if(all(tri) && { ULs <- sapply(Tlst, slot, "uplo")##  "U" or "L"
153      as(ret, "CsparseMatrix")                            all(ULs[1L] == ULs[-1L]) } ## all upper or all lower
154           ){ ## *triangular* result
155    
156            if(use.n) { ## return nsparseMatrix :
157                r <- new("ntTMatrix")
158            } else {
159                r <- new(paste0(if(use.l) "l" else "d", "tTMatrix"))
160                r@x <- unlist(lapply(Tlst, slot, "x"))
161            }
162            r@uplo <- ULs[1L]
163        }
164        else {
165            if(any(sym))
166                Tlst[sym] <- lapply(Tlst[sym], as, "generalMatrix")
167            if(use.n) { ## return nsparseMatrix :
168                r <- new("ngTMatrix")
169            } else {
170                r <- new(paste0(if(use.l) "l" else "d", "gTMatrix"))
171                r@x <- unlist(lapply(Tlst, slot, "x"))
172            }
173        }
174        r@Dim <- c(i_off[nl+1], j_off[nl + 1])
175        r@i <- unlist(lapply(1:nl, function(k) Tlst[[k]]@i + i_off[k]))
176        r@j <- unlist(lapply(1:nl, function(k) Tlst[[k]]@j + j_off[k]))
177        r
178    }
179    
180    bdiag <- function(...) {
181        if((nA <- nargs()) == 0) return(new("dgCMatrix"))
182        if(nA == 1 && !is.list(...))
183            return(as(..., "CsparseMatrix"))
184        alis <- if(nA == 1 && is.list(..1)) ..1 else list(...)
185        if(length(alis) == 1)
186            return(as(alis[[1]], "CsparseMatrix"))
187    
188        ## else : two or more arguments
189        as(.bdiag(alis), "CsparseMatrix")
190  }  }
191    
192  diag2T <- function(from) {  
193      i <- if(from@diag == "U") integer(0) else seq_len(from@Dim[1]) - 1:1  .diag2tT <- function(from, uplo = "U", kind = .M.kind(from)) {
194      new(paste(.M.kind(from), "tTMatrix", sep=''),      ## to triangular Tsparse
195        i <- if(from@diag == "U") integer(0) else seq_len(from@Dim[1]) - 1L
196        new(paste(kind, "tTMatrix", sep=''),
197          diag = from@diag, Dim = from@Dim, Dimnames = from@Dimnames,          diag = from@diag, Dim = from@Dim, Dimnames = from@Dimnames,
198            uplo = uplo,
199          x = from@x, # <- ok for diag = "U" and "N" (!)          x = from@x, # <- ok for diag = "U" and "N" (!)
200          i = i, j = i)          i = i, j = i)
201  }  }
202    
203  setAs("diagonalMatrix", "triangularMatrix", diag2T)  .diag2sT <- function(from, uplo = "U", kind = .M.kind(from)) {
204  setAs("diagonalMatrix", "sparseMatrix", diag2T)      ## to symmetric Tsparse
205        n <- from@Dim[1]
206        i <- seq_len(n) - 1L
207        new(paste(kind, "sTMatrix", sep=''),
208            Dim = from@Dim, Dimnames = from@Dimnames,
209            i = i, j = i, uplo = uplo,
210            x = if(from@diag == "N") from@x else ## "U"-diag
211            rep.int(switch(kind,
212                           "d" = 1.,
213                           "l" =,
214                           "n" = TRUE,
215                           ## otherwise
216                           stop("'", kind,"' kind not yet implemented")), n))
217    }
218    
219    ## diagonal -> triangular,  upper / lower depending on "partner":
220    diag2tT.u <- function(d, x, kind = .M.kind(d))
221        .diag2tT(d, uplo = if(is(x,"triangularMatrix")) x@uplo else "U", kind)
222    
223    ## diagonal -> sparse {triangular OR symmetric} (upper / lower) depending on "partner":
224    diag2Tsmart <- function(d, x, kind = .M.kind(d)) {
225        clx <- getClassDef(class(x))
226        if(extends(clx, "symmetricMatrix"))
227            .diag2sT(d, uplo = x@uplo, kind)
228        else
229            .diag2tT(d, uplo = if(extends(clx,"triangularMatrix")) x@uplo else "U", kind)
230    }
231    
232    
233    ## In order to evade method dispatch ambiguity warnings,
234    ## and because we can save a .M.kind() call, we use this explicit
235    ## "hack"  instead of signature  x = "diagonalMatrix" :
236    ##
237    ## ddi*:
238    diag2tT <- function(from) .diag2tT(from, "U", "d")
239    setAs("ddiMatrix", "triangularMatrix", diag2tT)
240    ##_no_longer_ setAs("ddiMatrix", "sparseMatrix", diag2tT)
241  ## needed too (otherwise <dense> -> Tsparse is taken):  ## needed too (otherwise <dense> -> Tsparse is taken):
242  setAs("diagonalMatrix", "TsparseMatrix", diag2T)  setAs("ddiMatrix", "TsparseMatrix", diag2tT)
243  ## is better than this:  setAs("ddiMatrix", "dsparseMatrix", diag2tT)
244  ## setAs("diagonalMatrix", "sparseMatrix",  setAs("ddiMatrix", "CsparseMatrix",
245  ##       function(from)        function(from) as(.diag2tT(from, "U", "d"), "CsparseMatrix"))
246  ##        as(from, if(is(from, "dMatrix")) "dgCMatrix" else "lgCMatrix"))  setAs("ddiMatrix", "symmetricMatrix",
247  setAs("diagonalMatrix", "CsparseMatrix",        function(from) .diag2sT(from, "U", "d"))
248        function(from) as(diag2T(from), "CsparseMatrix"))  ##
249    ## ldi*:
250    diag2tT <- function(from) .diag2tT(from, "U", "l")
251    setAs("ldiMatrix", "triangularMatrix", diag2tT)
252    ##_no_longer_ setAs("ldiMatrix", "sparseMatrix", diag2tT)
253    ## needed too (otherwise <dense> -> Tsparse is taken):
254    setAs("ldiMatrix", "TsparseMatrix", diag2tT)
255    setAs("ldiMatrix", "lsparseMatrix", diag2tT)
256    setAs("ldiMatrix", "CsparseMatrix",
257          function(from) as(.diag2tT(from, "U", "l"), "CsparseMatrix"))
258    setAs("ldiMatrix", "symmetricMatrix",
259          function(from) .diag2sT(from, "U", "l"))
260    
261  setAs("diagonalMatrix", "matrix",  
262    setAs("diagonalMatrix", "nMatrix",
263        function(from) {        function(from) {
264            n <- from@Dim[1]            n <- from@Dim[1]
265            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
266                                       } else from@x,            new("ntTMatrix", i = i, j = i, diag = from@diag,
267                 nrow = n, ncol = n)                Dim = from@Dim, Dimnames = from@Dimnames)
268        })        })
269    
270  setAs("diagonalMatrix", "generalMatrix", # prefer sparse:  setAs("diagonalMatrix", "nsparseMatrix", function(from) as(from, "nMatrix"))
       function(from) as(from, paste(.M.kind(from), "gCMatrix", sep='')))  
271    
272  .diag.x <- function(m) {  ## Cheap fast substitute for diag() which *does* preserve the mode of x :
273      if(m@diag == "U")  mkDiag <- function(x, n) {
274          rep.int(if(is.numeric(m@x)) 1. else TRUE,      y <- matrix(as0(mod=mode(x)), n,n)
275                  m@Dim[1])      if (n > 0) y[1L + 0:(n - 1L) * (n + 1L)] <- x
276      else m@x      y
 }  
   
 .diag.2N <- function(m) {  
     if(m@diag == "U") m@diag <- "N"  
     m  
277  }  }
278    
279  ## given the above, the following  4  coercions should be all unneeded;  setAs("diagonalMatrix", "matrix",
 ## we prefer triangular to general:  
 setAs("ddiMatrix", "dgTMatrix",  
280        function(from) {        function(from) {
281            .Deprecated("as(, \"sparseMatrix\")")            ## want "ldiMatrix" -> <logical> "matrix" :
282            n <- from@Dim[1]            mkDiag(if(from@diag == "U") as1(from@x) else from@x,
283            i <- seq_len(n) - 1:1                   n = from@Dim[1])
284            new("dgTMatrix", i = i, j = i, x = .diag.x(from),        })
               Dim = c(n,n), Dimnames = from@Dimnames) })  
285    
286  setAs("ddiMatrix", "dgCMatrix",  setMethod("as.vector", signature(x = "diagonalMatrix", mode="missing"),
287        function(from) as(as(from, "sparseMatrix"), "dgCMatrix"))            function(x, mode) {
288                  n <- x@Dim[1]
289                  mod.x <- mode(x@x)
290                  r <- vector(mod.x, length = n^2)
291                  if(n)
292                      r[1 + 0:(n - 1) * (n + 1)] <-
293                          if(x@diag == "U") as1(mod=mod.x) else x@x
294                  r
295              })
296    
297  setAs("ldiMatrix", "lgTMatrix",  setAs("diagonalMatrix", "generalMatrix", # prefer sparse:
298        function(from) {        function(from) as(as(from, "CsparseMatrix"), "generalMatrix"))
           .Deprecated("as(, \"sparseMatrix\")")  
           n <- from@Dim[1]  
           if(from@diag == "U") { # unit-diagonal  
               x <- rep.int(TRUE, n)  
               i <- seq_len(n) - 1:1  
           } else { # "normal"  
               nz <- nz.NA(from@x, na. = TRUE)  
               x <- from@x[nz]  
               i <- which(nz) - 1:1  
           }  
           new("lgTMatrix", i = i, j = i, x = x,  
               Dim = c(n,n), Dimnames = from@Dimnames) })  
299    
300  setAs("ldiMatrix", "lgCMatrix",  setAs("diagonalMatrix", "denseMatrix",
301        function(from) as(as(from, "lgTMatrix"), "lgCMatrix"))        function(from) as(as(from, "CsparseMatrix"), "denseMatrix"))
302    
303    .diag.x <- function(m) if(m@diag == "U") rep.int(as1(m@x), m@Dim[1]) else m@x
304    
305    .diag.2N <- function(m) {
306        if(m@diag == "U") m@diag <- "N"
307        m
308    }
309    
 if(FALSE) # now have faster  "ddense" -> "dge"  
310  setAs("ddiMatrix", "dgeMatrix",  setAs("ddiMatrix", "dgeMatrix",
311        function(from) as(as(from, "matrix"), "dgeMatrix"))        function(from) .Call(dup_mMatrix_as_dgeMatrix, from))
312    setAs("ddiMatrix", "ddenseMatrix",
313          function(from) as(as(from, "triangularMatrix"),"denseMatrix"))
314    setAs("ldiMatrix", "ldenseMatrix",
315          function(from) as(as(from, "triangularMatrix"),"denseMatrix"))
316    
317    
318  setAs("matrix", "diagonalMatrix",  setAs("matrix", "diagonalMatrix",
319        function(from) {        function(from) {
320            d <- dim(from)            d <- dim(from)
321            if(d[1] != (n <- d[2])) stop("non-square matrix")            if(d[1] != (n <- d[2])) stop("non-square matrix")
322            if(any(from[row(from) != col(from)] != 0))            if(any(from[row(from) != col(from)] != 0))
323                stop("has non-zero off-diagonal entries")                stop("matrix with non-zero off-diagonals cannot be coerced to diagonalMatrix")
324            x <- diag(from)            x <- diag(from)
325            if(is.logical(x)) {            if(is.logical(x)) {
326                cl <- "ldiMatrix"                cl <- "ldiMatrix"
327                uni <- all(x)                uni <- allTrue(x) ## uni := {is it unit-diagonal ?}
328            } else {            } else {
329                cl <- "ddiMatrix"                cl <- "ddiMatrix"
330                uni <- all(x == 1)                uni <- allTrue(x == 1)
331                storage.mode(x) <- "double"                storage.mode(x) <- "double"
332            } ## TODO: complex            } ## TODO: complex
333            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 344 
344            x <- diag(from)            x <- diag(from)
345            if(is.logical(x)) {            if(is.logical(x)) {
346                cl <- "ldiMatrix"                cl <- "ldiMatrix"
347                uni <- all(x)                uni <- allTrue(x)
348            } else {            } else {
349                cl <- "ddiMatrix"                cl <- "ddiMatrix"
350                uni <- all(x == 1)                uni <- allTrue(x == 1)
351                storage.mode(x) <- "double"                storage.mode(x) <- "double"
352            }            } ## TODO: complex
353            new(cl, Dim = c(n,n), diag = if(uni) "U" else "N",            new(cl, Dim = c(n,n), diag = if(uni) "U" else "N",
354                x = if(uni) x[FALSE] else x)                x = if(uni) x[FALSE] else x)
355        })        })
356    
357    
358  setMethod("diag", signature(x = "diagonalMatrix"),  setMethod("diag", signature(x = "diagonalMatrix"),
359            function(x = 1, nrow, ncol = n) .diag.x(x))            function(x = 1, nrow, ncol) .diag.x(x))
360    
361    subDiag <- function(x, i, j, ..., drop) {
362  subDiag <- function(x, i, j, drop) {      x <- as(x, "CsparseMatrix") ## << was "TsparseMatrix" (Csparse is faster now)
     x <- as(x, "sparseMatrix")  
363      x <- if(missing(i))      x <- if(missing(i))
364          x[, j, drop=drop]          x[, j, drop=drop]
365      else if(missing(j))      else if(missing(j))
366          x[i, , drop=drop]          if(nargs() == 4) x[i, , drop=drop] else x[i, drop=drop]
367      else      else
368          x[i,j, drop=drop]          x[i,j, drop=drop]
369      if(isDiagonal(x)) as(x, "diagonalMatrix") else x      if(isS4(x) && isDiagonal(x)) as(x, "diagonalMatrix") else x
370  }  }
371    
372  setMethod("[", signature(x = "diagonalMatrix", i = "index",  setMethod("[", signature(x = "diagonalMatrix", i = "index",
373                           j = "index", drop = "logical"), subDiag)                           j = "index", drop = "logical"), subDiag)
374  setMethod("[", signature(x = "diagonalMatrix", i = "index",  setMethod("[", signature(x = "diagonalMatrix", i = "index",
375                          j = "missing", drop = "logical"),                          j = "missing", drop = "logical"),
376            function(x, i, drop) subDiag(x, i=i, drop=drop))            function(x, i, j, ..., drop) {
377                  na <- nargs()
378                  Matrix.msg("diag[i,m,l] : nargs()=", na, .M.level = 2)
379                  if(na == 4)
380                       subDiag(x, i=i, , drop=drop)
381                  else subDiag(x, i=i,   drop=drop)
382              })
383  setMethod("[", signature(x = "diagonalMatrix", i = "missing",  setMethod("[", signature(x = "diagonalMatrix", i = "missing",
384                           j = "index", drop = "logical"),                           j = "index", drop = "logical"),
385            function(x, j, drop) subDiag(x, j=j, drop=drop))            function(x, i, j, ..., drop) subDiag(x, j=j, drop=drop))
386    
387  ## When you assign to a diagonalMatrix, the result should be  ## When you assign to a diagonalMatrix, the result should be
388  ## diagonal or sparse ---  ## diagonal or sparse ---
389  ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch  ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch
390  replDiag <- function(x, i, j, value) {  ## Only(?) current bug:  x[i] <- value  is wrong when  i is *vector*
391      x <- as(x, "sparseMatrix")  replDiag <- function(x, i, j, ..., value) {
392        x <- as(x, "TsparseMatrix")
393      if(missing(i))      if(missing(i))
394          x[, j] <- value          x[, j] <- value
395      else if(missing(j))      else if(missing(j)) { ##  x[i , ] <- v  *OR*   x[i] <- v
396            na <- nargs()
397    ##         message("diagnosing replDiag() -- nargs()= ", na)
398            if(na == 4)
399          x[i, ] <- value          x[i, ] <- value
400      else          else if(na == 3)
401                x[i] <- value
402            else stop("Internal bug: nargs()=",na,"; please report")
403        } else
404          x[i,j] <- value          x[i,j] <- value
405      if(isDiagonal(x)) as(x, "diagonalMatrix") else x      if(isDiagonal(x)) as(x, "diagonalMatrix") else x
406  }  }
407    
408  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index",
409                                  j = "index", value = "replValue"), replDiag)                                  j = "index", value = "replValue"), replDiag)
410    
411  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index",
412                                  j = "missing", value = "replValue"),                                  j = "missing", value = "replValue"),
413                   function(x, i, value) replDiag(x, i=i, value=value))                   function(x,i,j, ..., value) {
414                         ## message("before replDiag() -- nargs()= ", nargs())
415                         if(nargs() == 3)
416                             replDiag(x, i=i, value=value)
417                         else ## nargs() == 4 :
418                             replDiag(x, i=i, , value=value)
419                     })
420    
421    setReplaceMethod("[", signature(x = "diagonalMatrix",
422                                    i = "matrix", # 2-col.matrix
423                                    j = "missing", value = "replValue"),
424                     function(x,i,j, ..., value) {
425                         if(ncol(i) == 2) {
426                             if(all((ii <- i[,1]) == i[,2])) { # replace in diagonal only
427                                 if(x@diag == "U") {
428                                     one <- as1(x@x)
429                                     if(any(value != one | is.na(value))) {
430                                         x@diag <- "N"
431                                         x@x <- rep.int(one, x@Dim[1])
432                                     }
433                                 }
434                                 x@x[ii] <- value
435                                 x
436                             } else { ## no longer diagonal, but remain sparse:
437                                 x <- as(x, "TsparseMatrix")
438                                 x[i] <- value
439                                 x
440                             }
441                         }
442                         else { # behave as "base R": use as if vector
443                             x <- as(x, "matrix")
444                             x[i] <- value
445                             Matrix(x)
446                         }
447                     })
448    
449  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing",
450                                  j = "index", value = "replValue"),                                  j = "index", value = "replValue"),
451                   function(x, j, value) replDiag(x, j=j, value=value))                   function(x,i,j, ..., value) replDiag(x, j=j, value=value))
452    
453    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",
454                                    value = "sparseMatrix"),
455                     function (x, i, j, ..., value)
456                     callGeneric(x=x, , j=j, value = as(value, "sparseVector")))
457    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "missing",
458                                    value = "sparseMatrix"),
459                     function (x, i, j, ..., value)
460                     callGeneric(x=x, i=i, , value = as(value, "sparseVector")))
461    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "index",
462                                    value = "sparseMatrix"),
463                     function (x, i, j, ..., value)
464                     callGeneric(x=x, i=i, j=j, value = as(value, "sparseVector")))
465    
466    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",
467                                    value = "sparseVector"),
468                     replDiag)
469    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "missing",
470                                    value = "sparseVector"),
471                     replDiag)
472    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "index",
473                                    value = "sparseVector"),
474                     replDiag)
475    
476    
477  setMethod("t", signature(x = "diagonalMatrix"),  setMethod("t", signature(x = "diagonalMatrix"),
# Line 232  Line 482 
482  setMethod("isTriangular", signature(object = "diagonalMatrix"),  setMethod("isTriangular", signature(object = "diagonalMatrix"),
483            function(object) TRUE)            function(object) TRUE)
484  setMethod("isSymmetric", signature(object = "diagonalMatrix"),  setMethod("isSymmetric", signature(object = "diagonalMatrix"),
485            function(object) TRUE)            function(object, ...) TRUE)
486    
487  setMethod("chol", signature(x = "ddiMatrix"),# pivot = "ANY"  setMethod("symmpart", signature(x = "diagonalMatrix"), function(x) x)
488            function(x, pivot) {  setMethod("skewpart", signature(x = "diagonalMatrix"), setZero)
489                if(any(x@x < 0)) stop("chol() is undefined for diagonal matrix with negative entries")  
490    setMethod("chol", signature(x = "ddiMatrix"),
491              function(x, pivot, ...) {
492                  if(x@diag == "U") return(x)
493                  ## else
494                  if(any(x@x < 0))
495                      stop("chol() is undefined for diagonal matrix with negative entries")
496                x@x <- sqrt(x@x)                x@x <- sqrt(x@x)
497                x                x
498            })            })
499  ## chol(L) is L for logical diagonal:  ## chol(L) is L for logical diagonal:
500  setMethod("chol", signature(x = "ldiMatrix"), function(x, pivot) x)  setMethod("chol", signature(x = "ldiMatrix"), function(x, pivot, ...) x)
501    
502  setMethod("!", "ldiMatrix", function(e1) {  setMethod("determinant", signature(x = "diagonalMatrix", logarithm = "logical"),
503      if(e1@diag == "N")            function(x, logarithm, ...) mkDet(.diag.x(x), logarithm))
504          e1@x <- !e1@x  
505      else { ## "U"  setMethod("norm", signature(x = "diagonalMatrix", type = "character"),
506          e1@diag <- "N"            function(x, type, ...) {
507          e1@x <- rep.int(FALSE, e1@Dim[1])                if((n <- x@Dim[1]) == 0) return(0) # as for "sparseMatrix"
508                  type <- toupper(substr(type[1], 1, 1))
509                  isU <- (x@diag == "U") # unit-diagonal
510                  if(type == "F") sqrt(if(isU) n else sum(x@x^2))
511                  else { ## norm == "I","1","O","M" :
512                      if(isU) 1 else max(abs(x@x))
513      }      }
     e1  
514  })  })
515    
516    
517    
518  ## Basic Matrix Multiplication {many more to add}  ## Basic Matrix Multiplication {many more to add}
519  ##       ---------------------  ##       ---------------------
520  ## Note that "ldi" logical are treated as numeric  ## Note that "ldi" logical are treated as numeric
521  diagdiagprod <- function(x, y) {  diagdiagprod <- function(x, y) {
522      if(any(dim(x) != dim(y))) stop("non-matching dimensions")      dimCheck(x,y)
523      if(x@diag != "U") {      if(x@diag != "U") {
524          if(y@diag != "U") {          if(y@diag != "U") {
525              nx <- x@x * y@x              nx <- x@x * y@x
# Line 285  Line 547 
547    
548    
549  diagmatprod <- function(x, y) {  diagmatprod <- function(x, y) {
550        ## x is diagonalMatrix
551      dx <- dim(x)      dx <- dim(x)
552      dy <- dim(y)      dy <- dim(y)
553      if(dx[2] != dy[1]) stop("non-matching dimensions")      if(dx[2] != dy[1]) stop("non-matching dimensions")
     n <- dx[1]  
554      as(if(x@diag == "U") y else x@x * y, "Matrix")      as(if(x@diag == "U") y else x@x * y, "Matrix")
555  }  }
   
556  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),
557            diagmatprod)            diagmatprod)
558    ## sneaky .. :
559  formals(diagmatprod) <- alist(x=, y=NULL)  formals(diagmatprod) <- alist(x=, y=NULL)
560  setMethod("crossprod", signature(x = "diagonalMatrix", y = "matrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "matrix"),
561            diagmatprod)            diagmatprod)
 setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "matrix"),  
           diagmatprod)  
562    
563  diagdgeprod <- function(x, y) {  diagGeprod <- function(x, y) {
564      dx <- dim(x)      dx <- dim(x)
565      dy <- dim(y)      dy <- dim(y)
566      if(dx[2] != dy[1]) stop("non-matching dimensions")      if(dx[2] != dy[1]) stop("non-matching dimensions")
# Line 308  Line 568 
568          y@x <- x@x * y@x          y@x <- x@x * y@x
569      y      y
570  }  }
571  setMethod("%*%", signature(x = "diagonalMatrix", y = "dgeMatrix"),  setMethod("%*%", signature(x= "diagonalMatrix", y= "dgeMatrix"), diagGeprod)
572            diagdgeprod, valueClass = "dgeMatrix")  setMethod("%*%", signature(x= "diagonalMatrix", y= "lgeMatrix"), diagGeprod)
573  formals(diagdgeprod) <- alist(x=, y=NULL)  formals(diagGeprod) <- alist(x=, y=NULL)
574  setMethod("crossprod", signature(x = "diagonalMatrix", y = "dgeMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "dgeMatrix"),
575            diagdgeprod, valueClass = "dgeMatrix")            diagGeprod, valueClass = "dgeMatrix")
576    setMethod("crossprod", signature(x = "diagonalMatrix", y = "lgeMatrix"),
577              diagGeprod)
578    
579  setMethod("%*%", signature(x = "matrix", y = "diagonalMatrix"),  matdiagprod <- function(x, y) {
           function(x, y) {  
580                dx <- dim(x)                dx <- dim(x)
581                dy <- dim(y)                dy <- dim(y)
582                if(dx[2] != dy[1]) stop("non-matching dimensions")                if(dx[2] != dy[1]) stop("non-matching dimensions")
583                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]))
584            })  }
585    setMethod("%*%", signature(x = "matrix", y = "diagonalMatrix"),
586              matdiagprod)
587    formals(matdiagprod) <- alist(x=, y=NULL)
588    setMethod("tcrossprod", signature(x = "matrix", y = "diagonalMatrix"),
589              matdiagprod)
590    
591  setMethod("%*%", signature(x = "dgeMatrix", y = "diagonalMatrix"),  gediagprod <- function(x, y) {
           function(x, y) {  
592                dx <- dim(x)                dx <- dim(x)
593                dy <- dim(y)                dy <- dim(y)
594                if(dx[2] != dy[1]) stop("non-matching dimensions")                if(dx[2] != dy[1]) stop("non-matching dimensions")
595                if(y@diag == "N")                if(y@diag == "N")
596                    x@x <- x@x * rep(y@x, each = dx[1])                    x@x <- x@x * rep(y@x, each = dx[1])
597                x                x
598            })  }
599    setMethod("%*%", signature(x= "dgeMatrix", y= "diagonalMatrix"), gediagprod)
600    setMethod("%*%", signature(x= "lgeMatrix", y= "diagonalMatrix"), gediagprod)
601    formals(gediagprod) <- alist(x=, y=NULL)
602    setMethod("tcrossprod", signature(x = "dgeMatrix", y = "diagonalMatrix"),
603              gediagprod)
604    setMethod("tcrossprod", signature(x = "lgeMatrix", y = "diagonalMatrix"),
605              gediagprod)
606    
607  ## crossprod {more of these}  ## crossprod {more of these}
608    
609  ## tcrossprod --- all are not yet there: do the dense ones here:  ## tcrossprod --- all are not yet there: do the dense ones here:
610    
611    setMethod("%*%", signature(x = "diagonalMatrix", y = "denseMatrix"),
612              function(x, y) if(x@diag == "U") y else x %*% as(y, "generalMatrix"))
613    setMethod("%*%", signature(x = "denseMatrix", y = "diagonalMatrix"),
614              function(x, y) if(y@diag == "U") x else as(x, "generalMatrix") %*% y)
615    
616    
617  ## FIXME:  ## FIXME:
618  ## setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "denseMatrix"),  ## setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "denseMatrix"),
619  ##        function(x, y = NULL) {  ##        function(x, y = NULL) {
620  ##           })  ##           })
621    
622  ## setMethod("tcrossprod", signature(x = "denseMatrix", y = "diagonalMatrix"),  Cspdiagprod <- function(x, y) {
623  ##        function(x, y = NULL) {      dx <- dim(x <- .Call(Csparse_diagU2N, x))
624  ##           })      dy <- dim(y)
625        if(dx[2] != dy[1]) stop("non-matching dimensions")
626        if(y@diag == "N") {
627            if(!all(y@x[1L] == y@x[-1L]) && is(x, "symmetricMatrix"))
628                x <- as(x, "generalMatrix")
629            ind <- rep.int(seq_len(dx[2]), x@p[-1] - x@p[-dx[2]-1L])
630            x@x <- x@x * y@x[ind]
631        }
632        if(is(x, "compMatrix") && length(xf <- x@factors)) {
633            ## instead of dropping all factors, be smart about some
634            ## TODO ......
635            x@factors <- list()
636        }
637        x
638    }
639    
640    diagCspprod <- function(x, y) {
641        dx <- dim(x)
642        dy <- dim(y <- .Call(Csparse_diagU2N, y))
643        if(dx[2] != dy[1]) stop("non-matching dimensions")
644        if(x@diag == "N") {
645            if(!all(x@x[1L] == x@x[-1L]) && is(y, "symmetricMatrix"))
646                y <- as(y, "generalMatrix")
647            y@x <- y@x * x@x[y@i + 1L]
648        }
649        if(is(y, "compMatrix") && length(yf <- y@factors)) {
650            ## instead of dropping all factors, be smart about some
651            ## TODO
652            keep <- character()
653            if(iLU <- names(yf) == "LU") {
654                ## TODO keep <- "LU"
655            }
656            y@factors <- yf[keep]
657        }
658        y
659    }
660    
661    setMethod("crossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
662              function(x, y = NULL) diagCspprod(x, y))
663    
664  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
665            function(x, y = NULL) { x <- as(x, "sparseMatrix"); callGeneric() })            function(x, y = NULL) diagCspprod(x, as(y, "CsparseMatrix")))
666    
667    ## Prefer calling diagCspprod to Cspdiagprod if going to transpose anyway
668    ##  x'y == (y'x)'
669    setMethod("crossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
670              function(x, y = NULL) t(diagCspprod(y, x)))
671    
672  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
673            function(x, y = NULL) { y <- as(y, "sparseMatrix"); callGeneric() })            function(x, y = NULL) t(diagCspprod(y, as(x, "Csparsematrix"))))
674    
675    setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
676              function(x, y = NULL) diagCspprod(x, t(y)))
677    
678  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
679            function(x, y = NULL) { x <- as(x, "sparseMatrix"); callGeneric() })            function(x, y = NULL) diagCspprod(x, t(as(y, "CsparseMatrix"))))
680    
681    setMethod("tcrossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
682              function(x, y = NULL) Cspdiagprod(x, y))
683    
684  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
685            function(x, y = NULL) { y <- as(y, "sparseMatrix"); callGeneric() })            function(x, y = NULL) Cspdiagprod(as(x, "CsparseMatrix"), y))
686    
687    setMethod("%*%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
688              function(x, y) diagCspprod(x, y))
689    
 ## FIXME?: In theory, this can be done *FASTER*, in some cases, via tapply1()  
690  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
691            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)!  
692    
693  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
694            function(x, y) x %*% as(y, "sparseMatrix"))            function(x, y) Cspdiagprod(as(x, "CsparseMatrix"), y))
695    
696    setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
697              function(x, y) Cspdiagprod(x, y))
698    
699    ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*
700    ##       do indeed work by going through sparse (and *not* ddense)!
701    
702  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),
703            function(a, b, ...) {            function(a, b, ...) {
# Line 384  Line 707 
707            })            })
708    
709  solveDiag <- function(a, b, ...) {  solveDiag <- function(a, b, ...) {
710      if((n <- a@Dim[1]) != nrow(b))      if(a@Dim[1] != nrow(b))
711          stop("incompatible matrix dimensions")          stop("incompatible matrix dimensions")
712      ## trivially invert a 'in place' and multiply:      ## trivially invert a 'in place' and multiply:
713      a@x <- 1/ a@x      a@x <- 1/ a@x
# Line 396  Line 719 
719  setMethod("solve", signature(a = "diagonalMatrix", b = "Matrix"),  setMethod("solve", signature(a = "diagonalMatrix", b = "Matrix"),
720            solveDiag)            solveDiag)
721    
722    ## Schur()  ---> ./eigen.R
723    
724    
725    
726  ### ---------------- diagonal  o  sparse  -----------------------------  ###---------------- <Ops> (<Arith>, <Logic>, <Compare> ) ----------------------
   
727    
728  ## Use function for several signatures, in order to evade  ## Use function for several signatures, in order to evade
729  ## ambiguous dispatch for "ddi", since there's also Arith(ddense., ddense.)  ## ambiguous dispatch for "ddi", since there's also Arith(ddense., ddense.)
730  diagOdiag <- function(e1,e2) { # result should also be diagonal  diagOdiag <- function(e1,e2) {
731        ## result should also be diagonal _ if possible _
732      r <- callGeneric(.diag.x(e1), .diag.x(e2)) # error if not "compatible"      r <- callGeneric(.diag.x(e1), .diag.x(e2)) # error if not "compatible"
733        ## Check what happens with non-diagonals, i.e. (0 o 0), (FALSE o 0), ...:
734        r00 <- callGeneric(if(is.numeric(e1@x)) 0 else FALSE,
735                           if(is.numeric(e2@x)) 0 else FALSE)
736        if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* diagonal
737      if(is.numeric(r)) {      if(is.numeric(r)) {
738          if(is.numeric(e2@x)) {          if(is.numeric(e2@x)) {
739              e2@x <- r; return(.diag.2N(e2)) }              e2@x <- r; return(.diag.2N(e2)) }
# Line 419  Line 747 
747      e1@x <- r      e1@x <- r
748      .diag.2N(e1)      .diag.2N(e1)
749  }  }
750        else { ## result not diagonal, but at least symmetric:
751            ## e.g., m == m
752            isNum <- (is.numeric(r) || is.numeric(r00))
753            isLog <- (is.logical(r) || is.logical(r00))
754            Matrix.msg("exploding <diag> o <diag> into dense matrix", .M.level = 2)
755            d <- e1@Dim
756            n <- d[1]
757            stopifnot(length(r) == n)
758            xx <- as.vector(matrix(rbind(r, matrix(r00,n,n)), n,n))
759            newcl <-
760                paste(if(isNum) "d" else if(isLog) {
761                    if(!any(is.na(r)) && !any(is.na(r00))) "n" else "l"
762                } else stop("not yet implemented .. please report")
763                      ,
764                      "syMatrix", sep='')
765    
766            new(newcl, Dim = e1@Dim, Dimnames = e1@Dimnames, x = xx)
767        }
768    }
769    
770    ### This would be *the* way, but we get tons of "ambiguous method dispatch"
771    ## we use this hack instead of signature  x = "diagonalMatrix" :
772    diCls <- names(getClass("diagonalMatrix")@subclasses)
773    if(FALSE) {
774  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "diagonalMatrix"),  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "diagonalMatrix"),
775            diagOdiag)            diagOdiag)
776  ## These two are just for method disambiguation:  } else { ## These are just for method disambiguation:
777  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "diagonalMatrix"),      for(c1 in diCls)
778            diagOdiag)          for(c2 in diCls)
779  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "ddiMatrix"),              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)
780            diagOdiag)  }
781    
782    ## diagonal  o  triangular  |-->  triangular
783    ## diagonal  o  symmetric   |-->  symmetric
784    ##    {also when other is sparse: do these "here" --
785    ##     before conversion to sparse, since that loses "diagonality"}
786    
787  ## For almost everything else, diag* shall be treated "as sparse" :  ## For almost everything else, diag* shall be treated "as sparse" :
788  ## These are cheap implementations via coercion  ## These are cheap implementations via coercion
789    
790  ## for disambiguation  ## For disambiguation --- define this for "sparseMatrix" , then for "ANY";
791  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "sparseMatrix"),  ## and because we can save an .M.kind() call, we use this explicit
792            function(e1,e2) callGeneric(as(e1, "sparseMatrix"), e2))  ## "hack" for all diagonalMatrix *subclasses* instead of just "diagonalMatrix" :
793  setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "diagonalMatrix"),  ##
794            function(e1,e2) callGeneric(e1, as(e2, "sparseMatrix")))  ## ddi*:
795  ## in general:  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "sparseMatrix"),
796  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "ANY"),            function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2))
797            function(e1,e2) callGeneric(as(e1,"sparseMatrix"), e2))  setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "ddiMatrix"),
798  setMethod("Ops", signature(e1 = "ANY", e2 = "diagonalMatrix"),            function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "d")))
799            function(e1,e2) callGeneric(e1, as(e2,"sparseMatrix")))  ## ldi*
800    setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "sparseMatrix"),
801              function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2))
802    setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "ldiMatrix"),
803              function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "l")))
804    
805    ## Ops:  Arith  --> numeric : "dMatrix"
806    ##       Compare --> logical
807    ##       Logic   --> logical: "lMatrix"
808    
809    ##  other = "numeric" : stay diagonal if possible
810    ## ddi*: Arith: result numeric, potentially ddiMatrix
811    setMethod("Arith", signature(e1 = "ddiMatrix", e2 = "numeric"),
812              function(e1,e2) {
813                  n <- e1@Dim[1]; nsq <- n^2
814                  f0 <- callGeneric(0, e2)
815                  if(all(is0(f0))) { # remain diagonal
816                      L1 <- (le <- length(e2)) == 1L
817                      if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)
818                      if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {
819                          e1@diag <- "N"
820                          if(L1) r <- rep.int(r, n)
821                      } else
822                          r <- callGeneric(e1@x, e2)
823                      e1@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
824                      return(e1)
825                  }
826                  callGeneric(diag2tT.u(e1,e2, "d"), e2)
827              })
828    
829    setMethod("Arith", signature(e1 = "numeric", e2 = "ddiMatrix"),
830              function(e1,e2) {
831                  n <- e2@Dim[1]; nsq <- n^2
832                  f0 <- callGeneric(e1, 0)
833                  if(all(is0(f0))) { # remain diagonal
834                      L1 <- (le <- length(e1)) == 1L
835                      if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)
836                      if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {
837                          e2@diag <- "N"
838                          if(L1) r <- rep.int(r, n)
839                      } else
840                          r <- callGeneric(e1, e2@x)
841                      e2@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
842                      return(e2)
843                  }
844                  callGeneric(e1, diag2tT.u(e2,e1, "d"))
845              })
846    
847    ## ldi* Arith --> result numeric, potentially ddiMatrix
848    setMethod("Arith", signature(e1 = "ldiMatrix", e2 = "numeric"),
849              function(e1,e2) {
850                  n <- e1@Dim[1]; nsq <- n^2
851                  f0 <- callGeneric(0, e2)
852                  if(all(is0(f0))) { # remain diagonal
853                      L1 <- (le <- length(e2)) == 1L
854                      if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)
855                      if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {
856                          e1@diag <- "N"
857                          if(L1) r <- rep.int(r, n)
858                      } else
859                          r <- callGeneric(e1@x, e2)
860                      e1 <- copyClass(e1, "ddiMatrix", c("diag", "Dim", "Dimnames"))
861                      e1@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
862                      return(e1)
863                  }
864                  callGeneric(diag2tT.u(e1,e2, "d"), e2)
865              })
866    
867    setMethod("Arith", signature(e1 = "numeric", e2 = "ldiMatrix"),
868              function(e1,e2) {
869                  n <- e2@Dim[1]; nsq <- n^2
870                  f0 <- callGeneric(e1, 0)
871                  if(all(is0(f0))) { # remain diagonal
872                      L1 <- (le <- length(e1)) == 1L
873                      if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)
874                      if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {
875                          e2@diag <- "N"
876                          if(L1) r <- rep.int(r, n)
877                      } else
878                          r <- callGeneric(e1, e2@x)
879                      e2 <- copyClass(e2, "ddiMatrix", c("diag", "Dim", "Dimnames"))
880                      e2@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
881                      return(e2)
882                  }
883                  callGeneric(e1, diag2tT.u(e2,e1, "d"))
884              })
885    
886    ## ddi*: for "Ops" without Arith --> result logical, potentially ldi
887    setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "numeric"),
888              function(e1,e2) {
889                  n <- e1@Dim[1]; nsq <- n^2
890                  f0 <- callGeneric(0, e2)
891                  if(all(is0(f0))) { # remain diagonal
892                      L1 <- (le <- length(e2)) == 1L
893                      if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)
894                      if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {
895                          e1@diag <- "N"
896                          if(L1) r <- rep.int(r, n)
897                      } else
898                          r <- callGeneric(e1@x, e2)
899                      e1 <- copyClass(e1, "ldiMatrix", c("diag", "Dim", "Dimnames"))
900                      e1@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
901                      return(e1)
902                  }
903                  callGeneric(diag2tT.u(e1,e2, "l"), e2)
904              })
905    
906    setMethod("Ops", signature(e1 = "numeric", e2 = "ddiMatrix"),
907              function(e1,e2) {
908                  n <- e2@Dim[1]; nsq <- n^2
909                  f0 <- callGeneric(e1, 0)
910                  if(all(is0(f0))) { # remain diagonal
911                      L1 <- (le <- length(e1)) == 1L
912                      if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)
913                      if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {
914                          e2@diag <- "N"
915                          if(L1) r <- rep.int(r, n)
916                      } else
917                          r <- callGeneric(e1, e2@x)
918                      e2 <- copyClass(e2, "ldiMatrix", c("diag", "Dim", "Dimnames"))
919                      e2@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
920                      return(e2)
921                  }
922                  callGeneric(e1, diag2tT.u(e2,e1, "l"))
923              })
924    
925    ## ldi*: for "Ops" without Arith --> result logical, potentially ldi
926    setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "numeric"),
927              function(e1,e2) {
928                  n <- e1@Dim[1]; nsq <- n^2
929                  f0 <- callGeneric(FALSE, e2)
930                  if(all(is0(f0))) { # remain diagonal
931                      L1 <- (le <- length(e2)) == 1L
932                      if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)
933                      if(e1@diag == "U" && any((r <- callGeneric(TRUE, e2)) != 1)) {
934                          e1@diag <- "N"
935                          if(L1) r <- rep.int(r, n)
936                      } else
937                          r <- callGeneric(e1@x, e2)
938                      e1@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
939                      return(e1)
940                  }
941                  callGeneric(diag2tT.u(e1,e2, "l"), e2)
942              })
943    
944    setMethod("Ops", signature(e1 = "numeric", e2 = "ldiMatrix"),
945              function(e1,e2) {
946                  n <- e2@Dim[1]; nsq <- n^2
947                  f0 <- callGeneric(e1, FALSE)
948                  if(all(is0(f0))) { # remain diagonal
949                      L1 <- (le <- length(e1)) == 1L
950                      if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)
951                      if(e2@diag == "U" && any((r <- callGeneric(e1, TRUE)) != 1)) {
952                          e2@diag <- "N"
953                          if(L1) r <- rep.int(r, n)
954                      } else
955                          r <- callGeneric(e1, e2@x)
956                      e2@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
957                      return(e2)
958                  }
959                  callGeneric(e1, diag2tT.u(e2,e1, "l"))
960              })
961    
962    
963    
964    ## Not {"sparseMatrix", "numeric} :  {"denseMatrix", "matrix", ... }
965    for(other in c("ANY", "Matrix", "dMatrix")) {
966        ## ddi*:
967        setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other),
968                  function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2))
969        setMethod("Ops", signature(e1 = other, e2 = "ddiMatrix"),
970                  function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "d")))
971        ## ldi*:
972        setMethod("Ops", signature(e1 = "ldiMatrix", e2 = other),
973                  function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2))
974        setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),
975                  function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "l")))
976    }
977    
978    ## Direct subclasses of "denseMatrix": currently ddenseMatrix, ldense... :
979    if(FALSE)## too general, would contain  denseModelMatrix:
980    dense.subCl <- local({ dM.scl <- getClass("denseMatrix")@subclasses
981                           names(dM.scl)[sapply(dM.scl, slot, "distance") == 1] })
982    dense.subCl <- paste(c("d","l","n"), "denseMatrix", sep="")
983    for(DI in diCls) {
984        for(c2 in c(dense.subCl, "Matrix")) {
985            for(Fun in c("*", "^", "&")) {
986                setMethod(Fun, signature(e1 = DI, e2 = c2),
987                          function(e1,e2) callGeneric(e1, Diagonal(x = diag(e2))))
988                setMethod(Fun, signature(e1 = c2, e2 = DI),
989                          function(e1,e2) callGeneric(Diagonal(x = diag(e1)), e2))
990            }
991            ## NB: This arguably implicitly uses  0/0 :== 0  to keep diagonality
992            for(Fun in c("%%", "%/%", "/")) {
993                setMethod(Fun, signature(e1 = DI, e2 = c2),
994                          function(e1,e2) callGeneric(e1, Diagonal(x = diag(e2))))
995            }
996        }
997    }
998    
999    
1000    ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"
1001    ### ----------   the last 4: separately here
1002    for(cl in diCls) {
1003    setMethod("any", cl,
1004              function (x, ..., na.rm) {
1005                  if(any(x@Dim == 0)) FALSE
1006                  else if(x@diag == "U") TRUE else any(x@x, ..., na.rm = na.rm)
1007              })
1008    setMethod("all",  cl, function (x, ..., na.rm) {
1009        n <- x@Dim[1]
1010        if(n >= 2) FALSE
1011        else if(n == 0 || x@diag == "U") TRUE
1012        else all(x@x, ..., na.rm = na.rm)
1013    })
1014    setMethod("prod", cl, function (x, ..., na.rm) {
1015        n <- x@Dim[1]
1016        if(n >= 2) 0
1017        else if(n == 0 || x@diag == "U") 1
1018        else ## n == 1, diag = "N" :
1019            prod(x@x, ..., na.rm = na.rm)
1020    })
1021    
1022    setMethod("sum", cl,
1023              function(x, ..., na.rm) {
1024                  r <- sum(x@x, ..., na.rm = na.rm)# double or integer, correctly
1025                  if(x@diag == "U" && !is.na(r)) r + x@Dim[1] else r
1026              })
1027    }
1028    
1029    ## The remaining ones are  max, min, range :
1030    
1031    setMethod("Summary", "ddiMatrix",
1032              function(x, ..., na.rm) {
1033                  if(any(x@Dim == 0)) callGeneric(numeric(0), ..., na.rm=na.rm)
1034                  else if(x@diag == "U")
1035                      callGeneric(x@x, 0, 1, ..., na.rm=na.rm)
1036                  else callGeneric(x@x, 0, ..., na.rm=na.rm)
1037              })
1038    setMethod("Summary", "ldiMatrix",
1039              function(x, ..., na.rm) {
1040                  if(any(x@Dim == 0)) callGeneric(logical(0), ..., na.rm=na.rm)
1041                  else if(x@diag == "U")
1042                      callGeneric(x@x, FALSE, TRUE, ..., na.rm=na.rm)
1043                  else callGeneric(x@x, FALSE, ..., na.rm=na.rm)
1044              })
1045    
1046    
1047    
# Line 455  Line 1056 
1056      invisible(x)      invisible(x)
1057  }  }
1058    
1059    ## somewhat consistent with "print" for sparseMatrix :
1060    setMethod("print", signature(x = "diagonalMatrix"), prDiag)
1061    
1062  setMethod("show", signature(object = "diagonalMatrix"),  setMethod("show", signature(object = "diagonalMatrix"),
1063            function(object) {            function(object) {
1064                d <- dim(object)                d <- dim(object)
1065                cl <- class(object)                cl <- class(object)
1066                cat(sprintf('%d x %d diagonal matrix of class "%s"\n',                cat(sprintf('%d x %d diagonal matrix of class "%s"',
1067                            d[1], d[2], cl))                            d[1], d[2], cl))
1068                  if(d[1] < 50) {
1069                      cat("\n")
1070                prDiag(object)                prDiag(object)
1071                  } else {
1072                      cat(", with diagonal entries\n")
1073                      show(diag(object))
1074                      invisible(object)
1075                  }
1076            })            })

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

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