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

revision 2052, Wed Aug 15 13:33:19 2007 UTC revision 2209, Mon Jul 14 19:39:39 2008 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    ## Pkg 'spdep' had (relatively slow) versions of this as_dsCMatrix_I()
36    .symDiagonal <- function(n, x = rep.int(1,n), uplo = "U") {
37        stopifnot(n == (n. <- as.integer(n)), (n <- n.) >= 0)
38        if((lx <- length(x)) == 1) x <- rep.int(x, n)
39        else if(lx != n) stop("length(x) must be 1 or n")
40        cls <-
41            if(is.double(x)) "dsCMatrix"
42            else if(is.logical(x)) "lsCMatrix"
43            else { ## for now
44                storage.mode(x) <- "double"
45                "dsCMatrix"
46            }
47        new(cls, Dim = c(n,n), x = x, uplo = uplo,
48            i = if(n) 0:(n - 1L) else integer(0), p = 0:n)
49    }
50    
51  ### This is modified from a post of Bert Gunter to R-help on  1 Sep 2005.  ### This is modified from a post of Bert Gunter to R-help on  1 Sep 2005.
52  ### Bert's code built on a post by Andy Liaw who most probably was influenced  ### Bert's code built on a post by Andy Liaw who most probably was influenced
53  ### by earlier posts, notably one by Scott Chasalow on S-news, 16 Jan 2002  ### by earlier posts, notably one by Scott Chasalow on S-news, 16 Jan 2002
54  ### who posted his bdiag() function written in December 1995.  ### who posted his bdiag() function written in December 1995.
55    if(FALSE)##--- no longer used:
56  bdiag <- function(...) {  .bdiag <- function(lst) {
57      if(nargs() == 0) return(new("dgCMatrix"))      ### block-diagonal matrix [a dgTMatrix] from list of matrices
58      ## else :      stopifnot(is.list(lst), length(lst) >= 1)
59      mlist <- if (nargs() == 1) as.list(...) else list(...)      dims <- sapply(lst, dim, USE.NAMES=FALSE)
     dims <- sapply(mlist, dim)  
60      ## make sure we had all matrices:      ## make sure we had all matrices:
61      if(!(is.matrix(dims) && nrow(dims) == 2))      if(!(is.matrix(dims) && nrow(dims) == 2))
62          stop("some arguments are not matrices")          stop("some arguments are not matrices")
63      csdim <- rbind(rep.int(0L, 2),      csdim <- rbind(rep.int(0L, 2),
64                     apply(sapply(mlist, dim), 1, cumsum))                     apply(dims, 1, cumsum))
65      ret <- new("dgTMatrix", Dim = as.integer(csdim[nrow(csdim),]))      r <- new("dgTMatrix")
66        r@Dim <- as.integer(csdim[nrow(csdim),])
67      add1 <- matrix(1:0, 2,2)      add1 <- matrix(1:0, 2,2)
68      for(i in seq_along(mlist)) {      for(i in seq_along(lst)) {
69          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])
70          if(is.null(dim(indx))) ## non-square matrix          if(is.null(dim(indx))) ## non-square matrix
71              ret[indx[[1]],indx[[2]]] <- mlist[[i]]              r[indx[[1]],indx[[2]]] <- lst[[i]]
72          else ## square matrix          else ## square matrix
73              ret[indx[,1],indx[,2]] <- mlist[[i]]              r[indx[,1], indx[,2]] <- lst[[i]]
74        }
75        r
76    }
77    ## expand(<mer>) needed something like bdiag() for lower-triangular
78    ## (Tsparse) Matrices; hence Doug Bates provided a much more efficient
79    ##  implementation for those; now extended and generalized:
80    .bdiag <- function(lst) {
81        ## block-diagonal matrix [a dgTMatrix] from list of matrices
82        stopifnot(is.list(lst), (nl <- length(lst)) >= 1)
83    
84        Tlst <- lapply(lapply(lst, Matrix:::as_Csp2), # includes "diagU2N"
85                       as, "TsparseMatrix")
86        if(nl == 1) return(Tlst[[1]])
87        ## else
88        i_off <- c(0L, cumsum(sapply(Tlst, nrow)))
89        j_off <- c(0L, cumsum(sapply(Tlst, ncol)))
90    
91        clss <- sapply(Tlst, class)
92        knds <- substr(clss, 2, 2)
93        sym  <- knds == "s" # symmetric ones
94        tri  <- knds == "t" # triangular ones
95        use.n <- any(is.n <- substr(clss,1,1) == "n")
96        if(use.n && !(use.n <- all(is.n)))
97            Tlst[is.n] <- lapply(Tlst[is.n], as, "lMatrix")
98        if(all(sym)) { ## result should be *symmetric*
99            uplos <- sapply(Tlst, slot, "uplo") ## either "U" or "L"
100            tLU <- table(uplos)# of length 1 or 2 ..
101            if(length(tLU) == 1) { ## all "U" or all "L"
102                useU <- uplos[1] == "U"
103            } else { ## length(tLU) == 2, counting "L" and "U"
104                useU <- diff(tLU) >= 0
105                if(useU && (hasL <- tLU[1] > 0))
106                    Tlst[hasL] <- lapply(Tlst[hasL], t)
107                else if(!useU && (hasU <- tLU[2] > 0))
108                    Tlst[hasU] <- lapply(Tlst[hasU], t)
109            }
110            if(use.n) { ## return nsparseMatrix :
111                r <- new("nsTMatrix")
112            } else {
113                r <- new("dsTMatrix")
114                r@x <- unlist(lapply(Tlst, slot, "x"))
115            }
116            r@uplo <- if(useU) "U" else "L"
117        }
118        else if(all(tri) && { ULs <- sapply(Tlst, slot, "uplo")##  "U" or "L"
119                              all(ULs[1L] == ULs[-1L]) } ## all upper or all lower
120           ){ ## *triangular* result
121    
122            if(use.n) { ## return nsparseMatrix :
123                r <- new("ntTMatrix")
124            } else {
125                r <- new("dtTMatrix")
126                r@x <- unlist(lapply(Tlst, slot, "x"))
127            }
128            r@uplo <- ULs[1L]
129        }
130        else {
131            if(any(sym))
132                Tlst[sym] <- lapply(Tlst[sym], as, "generalMatrix")
133            if(use.n) { ## return nsparseMatrix :
134                r <- new("ngTMatrix")
135            } else {
136                r <- new("dgTMatrix")
137                r@x <- unlist(lapply(Tlst, slot, "x"))
138            }
139      }      }
140      ## slightly debatable if we really should return Csparse.. :      r@Dim <- c(i_off[nl+1], j_off[nl + 1])
141      as(ret, "CsparseMatrix")      r@i <- unlist(lapply(1:nl, function(k) Tlst[[k]]@i + i_off[k]))
142        r@j <- unlist(lapply(1:nl, function(k) Tlst[[k]]@j + j_off[k]))
143        r
144  }  }
145    
146  diag2tT <- function(from) {  bdiag <- function(...) {
147        if((nA <- nargs()) == 0) return(new("dgCMatrix"))
148        if(nA == 1 && !is.list(...))
149            return(as(..., "CsparseMatrix"))
150        alis <- if(nA == 1 && is.list(..1)) ..1 else list(...)
151        if(length(alis) == 1)
152            return(as(alis[[1]], "CsparseMatrix"))
153    
154        ## else : two or more arguments
155        as(.bdiag(alis), "CsparseMatrix")
156    }
157    
158    
159    .diag2tT <- function(from, uplo = "U", kind = .M.kind(from)) {
160        ## to triangular Tsparse
161      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
162      new(paste(.M.kind(from), "tTMatrix", sep=''),      new(paste(kind, "tTMatrix", sep=''),
163          diag = from@diag, Dim = from@Dim, Dimnames = from@Dimnames,          diag = from@diag, Dim = from@Dim, Dimnames = from@Dimnames,
164            uplo = uplo,
165          x = from@x, # <- ok for diag = "U" and "N" (!)          x = from@x, # <- ok for diag = "U" and "N" (!)
166          i = i, j = i)          i = i, j = i)
167  }  }
168    
169  diag2sT <- function(from) { # to symmetric Tsparse  .diag2sT <- function(from, uplo = "U", kind = .M.kind(from)) {
170      i <- if(from@diag == "U") integer(0) else seq_len(from@Dim[1]) - 1L      ## to symmetric Tsparse
171      new(paste(.M.kind(from), "sTMatrix", sep=''),      n <- from@Dim[1]
172        i <- seq_len(n) - 1L
173        new(paste(kind, "sTMatrix", sep=''),
174          Dim = from@Dim, Dimnames = from@Dimnames,          Dim = from@Dim, Dimnames = from@Dimnames,
175          x = from@x, i = i, j = i)          i = i, j = i, uplo = uplo,
176            x = if(from@diag == "N") from@x else ## "U"-diag
177            rep.int(switch(kind,
178                           "d" = 1.,
179                           "l" =,
180                           "n" = TRUE,
181                           ## otherwise
182                           stop("'", kind,"' kind not yet implemented")), n))
183    }
184    
185    ## diagonal -> triangular,  upper / lower depending on "partner":
186    diag2tT.u <- function(d, x, kind = .M.kind(d))
187        .diag2tT(d, uplo = if(is(x,"triangularMatrix")) x@uplo else "U", kind)
188    
189    ## diagonal -> sparse {triangular OR symmetric} (upper / lower) depending on "partner":
190    diag2Tsmart <- function(d, x, kind = .M.kind(d)) {
191        clx <- getClassDef(class(x))
192        if(extends(clx, "symmetricMatrix"))
193            .diag2sT(d, uplo = x@uplo, kind)
194        else
195            .diag2tT(d, uplo = if(extends(clx,"triangularMatrix")) x@uplo else "U", kind)
196  }  }
197    
198  setAs("diagonalMatrix", "triangularMatrix", diag2tT)  
199  setAs("diagonalMatrix", "sparseMatrix", diag2tT)  ## In order to evade method dispatch ambiguity warnings,
200    ## and because we can save a .M.kind() call, we use this explicit
201    ## "hack"  instead of signature  x = "diagonalMatrix" :
202    ##
203    ## ddi*:
204    diag2tT <- function(from) .diag2tT(from, "U", "d")
205    setAs("ddiMatrix", "triangularMatrix", diag2tT)
206    setAs("ddiMatrix", "sparseMatrix", diag2tT)
207  ## needed too (otherwise <dense> -> Tsparse is taken):  ## needed too (otherwise <dense> -> Tsparse is taken):
208  setAs("diagonalMatrix", "TsparseMatrix", diag2tT)  setAs("ddiMatrix", "TsparseMatrix", diag2tT)
209  ## is better than this:  setAs("ddiMatrix", "CsparseMatrix",
210  ## setAs("diagonalMatrix", "sparseMatrix",        function(from) as(.diag2tT(from, "U", "d"), "CsparseMatrix"))
211  ##       function(from)  setAs("ddiMatrix", "symmetricMatrix",
212  ##        as(from, if(is(from, "dMatrix")) "dgCMatrix" else "lgCMatrix"))        function(from) .diag2sT(from, "U", "d"))
213  setAs("diagonalMatrix", "CsparseMatrix",  ##
214        function(from) as(diag2tT(from), "CsparseMatrix"))  ## ldi*:
215    diag2tT <- function(from) .diag2tT(from, "U", "l")
216    setAs("ldiMatrix", "triangularMatrix", diag2tT)
217    setAs("ldiMatrix", "sparseMatrix", diag2tT)
218    ## needed too (otherwise <dense> -> Tsparse is taken):
219    setAs("ldiMatrix", "TsparseMatrix", diag2tT)
220    setAs("ldiMatrix", "CsparseMatrix",
221          function(from) as(.diag2tT(from, "U", "l"), "CsparseMatrix"))
222    setAs("ldiMatrix", "symmetricMatrix",
223          function(from) .diag2sT(from, "U", "l"))
224    
 setAs("diagonalMatrix", "symmetricMatrix", diag2sT)  
225    
226  setAs("diagonalMatrix", "matrix",  setAs("diagonalMatrix", "nMatrix",
227        function(from) {        function(from) {
228            n <- from@Dim[1]            n <- from@Dim[1]
229            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
230                                       } else from@x,            new("ntTMatrix", i = i, j = i, diag = from@diag,
231                 nrow = n, ncol = n)                Dim = from@Dim, Dimnames = from@Dimnames)
232        })        })
233    
 setAs("diagonalMatrix", "generalMatrix", # prefer sparse:  
       function(from) as(from, paste(.M.kind(from), "gCMatrix", sep='')))  
   
 .diag.x <- function(m) {  
     if(m@diag == "U")  
         rep.int(if(is.numeric(m@x)) 1. else TRUE,  
                 m@Dim[1])  
     else m@x  
 }  
234    
235  .diag.2N <- function(m) {  ## Cheap fast substitute for diag() which *does* preserve the mode of x :
236      if(m@diag == "U") m@diag <- "N"  mkDiag <- function(x, n) {
237      m      y <- matrix(as0(mod=mode(x)), n,n)
238        if (n > 0) y[1L + 0:(n - 1L) * (n + 1L)] <- x
239        y
240  }  }
241    
242  ## given the above, the following  4  coercions should be all unneeded;  setAs("diagonalMatrix", "matrix",
 ## we prefer triangular to general:  
 setAs("ddiMatrix", "dgTMatrix",  
243        function(from) {        function(from) {
244            .Deprecated("as(, \"sparseMatrix\")")            ## want "ldiMatrix" -> <logical> "matrix" :
245            n <- from@Dim[1]            mkDiag(if(from@diag == "U") as1(from@x) else from@x,
246            i <- seq_len(n) - 1L                   n = from@Dim[1])
247            new("dgTMatrix", i = i, j = i, x = .diag.x(from),        })
               Dim = c(n,n), Dimnames = from@Dimnames) })  
248    
249  setAs("ddiMatrix", "dgCMatrix",  setMethod("as.vector", signature(x = "diagonalMatrix", mode="missing"),
250        function(from) as(as(from, "sparseMatrix"), "dgCMatrix"))            function(x, mode) {
251                  n <- x@Dim[1]
252                  mod.x <- mode(x@x)
253                  r <- vector(mod.x, length = n^2)
254                  if(n)
255                      r[1 + 0:(n - 1) * (n + 1)] <-
256                          if(x@diag == "U") as1(mod=mod.x) else x@x
257                  r
258              })
259    
260  setAs("ldiMatrix", "lgTMatrix",  setAs("diagonalMatrix", "generalMatrix", # prefer sparse:
261        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) - 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) })  
262    
263  setAs("ldiMatrix", "lgCMatrix",  .diag.x <- function(m) if(m@diag == "U") rep.int(as1(m@x), m@Dim[1]) else m@x
       function(from) as(as(from, "lgTMatrix"), "lgCMatrix"))  
264    
265    .diag.2N <- function(m) {
266        if(m@diag == "U") m@diag <- "N"
267        m
268    }
269    
 if(FALSE) # now have faster  "ddense" -> "dge"  
270  setAs("ddiMatrix", "dgeMatrix",  setAs("ddiMatrix", "dgeMatrix",
271        function(from) as(as(from, "matrix"), "dgeMatrix"))        function(from) .Call(dup_mMatrix_as_dgeMatrix, from))
272    setAs("ddiMatrix", "ddenseMatrix",
273          function(from) as(as(from, "triangularMatrix"),"denseMatrix"))
274    setAs("ldiMatrix", "ldenseMatrix",
275          function(from) as(as(from, "triangularMatrix"),"denseMatrix"))
276    
277    
278  setAs("matrix", "diagonalMatrix",  setAs("matrix", "diagonalMatrix",
279        function(from) {        function(from) {
# Line 188  Line 318 
318  setMethod("diag", signature(x = "diagonalMatrix"),  setMethod("diag", signature(x = "diagonalMatrix"),
319            function(x = 1, nrow, ncol) .diag.x(x))            function(x = 1, nrow, ncol) .diag.x(x))
320    
321    subDiag <- function(x, i, j, ..., drop) {
 subDiag <- function(x, i, j, drop) {  
322      x <- as(x, "sparseMatrix")      x <- as(x, "sparseMatrix")
323      x <- if(missing(i))      x <- if(missing(i))
324          x[, j, drop=drop]          x[, j, drop=drop]
# Line 197  Line 326 
326          x[i, , drop=drop]          x[i, , drop=drop]
327      else      else
328          x[i,j, drop=drop]          x[i,j, drop=drop]
329      if(isDiagonal(x)) as(x, "diagonalMatrix") else x      if(isS4(x) && isDiagonal(x)) as(x, "diagonalMatrix") else x
330  }  }
331    
332  setMethod("[", signature(x = "diagonalMatrix", i = "index",  setMethod("[", signature(x = "diagonalMatrix", i = "index",
333                           j = "index", drop = "logical"), subDiag)                           j = "index", drop = "logical"), subDiag)
334  setMethod("[", signature(x = "diagonalMatrix", i = "index",  setMethod("[", signature(x = "diagonalMatrix", i = "index",
335                          j = "missing", drop = "logical"),                          j = "missing", drop = "logical"),
336            function(x, i, drop) subDiag(x, i=i, drop=drop))            function(x, i, j, ..., drop) subDiag(x, i=i, drop=drop))
337  setMethod("[", signature(x = "diagonalMatrix", i = "missing",  setMethod("[", signature(x = "diagonalMatrix", i = "missing",
338                           j = "index", drop = "logical"),                           j = "index", drop = "logical"),
339            function(x, j, drop) subDiag(x, j=j, drop=drop))            function(x, i, j, ..., drop) subDiag(x, j=j, drop=drop))
340    
341  ## When you assign to a diagonalMatrix, the result should be  ## When you assign to a diagonalMatrix, the result should be
342  ## diagonal or sparse ---  ## diagonal or sparse ---
343  ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch  ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch
344  replDiag <- function(x, i, j, value) {  ## Only(?) current bug:  x[i] <- value  is wrong when  i is *vector*
345    replDiag <- function(x, i, j, ..., value) {
346      x <- as(x, "sparseMatrix")      x <- as(x, "sparseMatrix")
347      if(missing(i))      if(missing(i))
348          x[, j] <- value          x[, j] <- value
349      else if(missing(j))      else if(missing(j)) { ##  x[i , ] <- v  *OR*   x[i] <- v
350            na <- nargs()
351    ##         message("diagnosing replDiag() -- nargs()= ", na)
352            if(na == 4)
353          x[i, ] <- value          x[i, ] <- value
354      else          else if(na == 3)
355                x[i] <- value
356            else stop("Internal bug: nargs()=",na,"; please report")
357        } else
358          x[i,j] <- value          x[i,j] <- value
359      if(isDiagonal(x)) as(x, "diagonalMatrix") else x      if(isDiagonal(x)) as(x, "diagonalMatrix") else x
360  }  }
361    
362  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index",
363                                  j = "index", value = "replValue"), replDiag)                                  j = "index", value = "replValue"), replDiag)
364    
365  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index",
366                                  j = "missing", value = "replValue"),                                  j = "missing", value = "replValue"),
367                   function(x, i, value) replDiag(x, i=i, value=value))                   function(x,i,j, ..., value) {
368                         ## message("before replDiag() -- nargs()= ", nargs())
369                         if(nargs() == 3)
370                             replDiag(x, i=i, value=value)
371                         else ## nargs() == 4 :
372                             replDiag(x, i=i, , value=value)
373                     })
374    
375    setReplaceMethod("[", signature(x = "diagonalMatrix",
376                                    i = "matrix", # 2-col.matrix
377                                    j = "missing", value = "replValue"),
378                     function(x,i,j, ..., value) {
379                         if(ncol(i) == 2) {
380                             if(all((ii <- i[,1]) == i[,2])) { # replace in diagonal only
381                                 if(x@diag == "U") {
382                                     one <- as1(x@x)
383                                     if(any(value != one | is.na(value))) {
384                                         x@diag <- "N"
385                                         x@x <- rep.int(one, x@Dim[1])
386                                     }
387                                 }
388                                 x@x[ii] <- value
389                                 x
390                             } else { ## no longer diagonal, but remain sparse:
391                                 x <- as(x, "sparseMatrix")
392                                 x[i] <- value
393                                 x
394                             }
395                         }
396                         else { # behave as "base R": use as if vector
397                             x <- as(x, "matrix")
398                             x[i] <- value
399                             Matrix(x)
400                         }
401                     })
402    
403  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing",
404                                  j = "index", value = "replValue"),                                  j = "index", value = "replValue"),
405                   function(x, j, value) replDiag(x, j=j, value=value))                   function(x,i,j, ..., value) replDiag(x, j=j, value=value))
406    
407    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",
408                                    value = "scarceMatrix"),
409                     function (x, i, j, ..., value)
410                     callGeneric(x=x, , j=j, value = as(value, "sparseVector")))
411    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "missing",
412                                    value = "scarceMatrix"),
413                     function (x, i, j, ..., value)
414                     callGeneric(x=x, i=i, , value = as(value, "sparseVector")))
415    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "index",
416                                    value = "scarceMatrix"),
417                     function (x, i, j, ..., value)
418                     callGeneric(x=x, i=i, j=j, value = as(value, "sparseVector")))
419    
420    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",
421                                    value = "sparseVector"),
422                     replDiag)
423    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "missing",
424                                    value = "sparseVector"),
425                     replDiag)
426    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "index",
427                                    value = "sparseVector"),
428                     replDiag)
429    
430    
431  setMethod("t", signature(x = "diagonalMatrix"),  setMethod("t", signature(x = "diagonalMatrix"),
# Line 241  Line 436 
436  setMethod("isTriangular", signature(object = "diagonalMatrix"),  setMethod("isTriangular", signature(object = "diagonalMatrix"),
437            function(object) TRUE)            function(object) TRUE)
438  setMethod("isSymmetric", signature(object = "diagonalMatrix"),  setMethod("isSymmetric", signature(object = "diagonalMatrix"),
439            function(object) TRUE)            function(object, ...) TRUE)
440    
441  setMethod("chol", signature(x = "ddiMatrix"),# pivot = "ANY"  setMethod("symmpart", signature(x = "diagonalMatrix"), function(x) x)
442            function(x, pivot) {  setMethod("skewpart", signature(x = "diagonalMatrix"), setZero)
443                if(any(x@x < 0)) stop("chol() is undefined for diagonal matrix with negative entries")  
444    setMethod("chol", signature(x = "ddiMatrix"),
445              function(x, pivot, ...) {
446                  if(x@diag == "U") return(x)
447                  ## else
448                  if(any(x@x < 0))
449                      stop("chol() is undefined for diagonal matrix with negative entries")
450                x@x <- sqrt(x@x)                x@x <- sqrt(x@x)
451                x                x
452            })            })
453  ## chol(L) is L for logical diagonal:  ## chol(L) is L for logical diagonal:
454  setMethod("chol", signature(x = "ldiMatrix"), function(x, pivot) x)  setMethod("chol", signature(x = "ldiMatrix"), function(x, pivot, ...) x)
455    
456    setMethod("determinant", signature(x = "diagonalMatrix", logarithm = "logical"),
457              function(x, logarithm, ...) mkDet(.diag.x(x), logarithm))
458    
459    setMethod("norm", signature(x = "diagonalMatrix", type = "character"),
460              function(x, type, ...) {
461                  if((n <- x@Dim[1]) == 0) return(0) # as for "sparseMatrix"
462                  type <- toupper(substr(type[1], 1, 1))
463                  isU <- (x@diag == "U") # unit-diagonal
464                  if(type == "F") sqrt(if(isU) n else sum(x@x^2))
465                  else { ## norm == "I","1","O","M" :
466                      if(isU) 1 else max(abs(x@x))
467                  }
468              })
469    
470    
471    
472  ## Basic Matrix Multiplication {many more to add}  ## Basic Matrix Multiplication {many more to add}
473  ##       ---------------------  ##       ---------------------
474  ## Note that "ldi" logical are treated as numeric  ## Note that "ldi" logical are treated as numeric
475  diagdiagprod <- function(x, y) {  diagdiagprod <- function(x, y) {
476      if(any(dim(x) != dim(y))) stop("non-matching dimensions")      n <- dimCheck(x,y)[1]
477      if(x@diag != "U") {      if(x@diag != "U") {
478          if(y@diag != "U") {          if(y@diag != "U") {
479              nx <- x@x * y@x              nx <- x@x * y@x
# Line 284  Line 501 
501    
502    
503  diagmatprod <- function(x, y) {  diagmatprod <- function(x, y) {
504        ## x is diagonalMatrix
505      dx <- dim(x)      dx <- dim(x)
506      dy <- dim(y)      dy <- dim(y)
507      if(dx[2] != dy[1]) stop("non-matching dimensions")      if(dx[2] != dy[1]) stop("non-matching dimensions")
508      n <- dx[1]      n <- dx[1]
509      as(if(x@diag == "U") y else x@x * y, "Matrix")      as(if(x@diag == "U") y else x@x * y, "Matrix")
510  }  }
   
511  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),
512            diagmatprod)            diagmatprod)
513    ## sneaky .. :
514  formals(diagmatprod) <- alist(x=, y=NULL)  formals(diagmatprod) <- alist(x=, y=NULL)
515  setMethod("crossprod", signature(x = "diagonalMatrix", y = "matrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "matrix"),
516            diagmatprod)            diagmatprod)
 setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "matrix"),  
           diagmatprod)  
517    
518  diagdgeprod <- function(x, y) {  diagGeprod <- function(x, y) {
519      dx <- dim(x)      dx <- dim(x)
520      dy <- dim(y)      dy <- dim(y)
521      if(dx[2] != dy[1]) stop("non-matching dimensions")      if(dx[2] != dy[1]) stop("non-matching dimensions")
# Line 307  Line 523 
523          y@x <- x@x * y@x          y@x <- x@x * y@x
524      y      y
525  }  }
526  setMethod("%*%", signature(x = "diagonalMatrix", y = "dgeMatrix"),  setMethod("%*%", signature(x= "diagonalMatrix", y= "dgeMatrix"), diagGeprod)
527            diagdgeprod, valueClass = "dgeMatrix")  setMethod("%*%", signature(x= "diagonalMatrix", y= "lgeMatrix"), diagGeprod)
528  formals(diagdgeprod) <- alist(x=, y=NULL)  formals(diagGeprod) <- alist(x=, y=NULL)
529  setMethod("crossprod", signature(x = "diagonalMatrix", y = "dgeMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "dgeMatrix"),
530            diagdgeprod, valueClass = "dgeMatrix")            diagGeprod, valueClass = "dgeMatrix")
531    setMethod("crossprod", signature(x = "diagonalMatrix", y = "lgeMatrix"),
532              diagGeprod)
533    
534  setMethod("%*%", signature(x = "matrix", y = "diagonalMatrix"),  matdiagprod <- function(x, y) {
           function(x, y) {  
535                dx <- dim(x)                dx <- dim(x)
536                dy <- dim(y)                dy <- dim(y)
537                if(dx[2] != dy[1]) stop("non-matching dimensions")                if(dx[2] != dy[1]) stop("non-matching dimensions")
538                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]))
539            })  }
540    setMethod("%*%", signature(x = "matrix", y = "diagonalMatrix"),
541              matdiagprod)
542    formals(matdiagprod) <- alist(x=, y=NULL)
543    setMethod("tcrossprod", signature(x = "matrix", y = "diagonalMatrix"),
544              matdiagprod)
545    
546  setMethod("%*%", signature(x = "dgeMatrix", y = "diagonalMatrix"),  gediagprod <- function(x, y) {
           function(x, y) {  
547                dx <- dim(x)                dx <- dim(x)
548                dy <- dim(y)                dy <- dim(y)
549                if(dx[2] != dy[1]) stop("non-matching dimensions")                if(dx[2] != dy[1]) stop("non-matching dimensions")
550                if(y@diag == "N")                if(y@diag == "N")
551                    x@x <- x@x * rep(y@x, each = dx[1])                    x@x <- x@x * rep(y@x, each = dx[1])
552                x                x
553            })  }
554    setMethod("%*%", signature(x= "dgeMatrix", y= "diagonalMatrix"), gediagprod)
555    setMethod("%*%", signature(x= "lgeMatrix", y= "diagonalMatrix"), gediagprod)
556    formals(gediagprod) <- alist(x=, y=NULL)
557    setMethod("tcrossprod", signature(x = "dgeMatrix", y = "diagonalMatrix"),
558              gediagprod)
559    setMethod("tcrossprod", signature(x = "lgeMatrix", y = "diagonalMatrix"),
560              gediagprod)
561    
562  ## crossprod {more of these}  ## crossprod {more of these}
563    
564  ## tcrossprod --- all are not yet there: do the dense ones here:  ## tcrossprod --- all are not yet there: do the dense ones here:
565    
566    setMethod("%*%", signature(x = "diagonalMatrix", y = "denseMatrix"),
567              function(x, y) if(x@diag == "U") y else x %*% as(y, "generalMatrix"))
568    setMethod("%*%", signature(x = "denseMatrix", y = "diagonalMatrix"),
569              function(x, y) if(y@diag == "U") x else as(x, "generalMatrix") %*% y)
570    
571    
572  ## FIXME:  ## FIXME:
573  ## setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "denseMatrix"),  ## setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "denseMatrix"),
574  ##        function(x, y = NULL) {  ##        function(x, y = NULL) {
# Line 345  Line 579 
579  ##           })  ##           })
580    
581  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
582            function(x, y = NULL) { x <- as(x, "sparseMatrix"); callGeneric() })            function(x, y = NULL) crossprod(as(x, "sparseMatrix"), y))
583    
584  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
585            function(x, y = NULL) { y <- as(y, "sparseMatrix"); callGeneric() })            function(x, y = NULL) crossprod(x, as(y, "sparseMatrix")))
586    
587  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
588            function(x, y = NULL) { x <- as(x, "sparseMatrix"); callGeneric() })            function(x, y = NULL) tcrossprod(as(x, "sparseMatrix"), y))
589    
590  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
591            function(x, y = NULL) { y <- as(y, "sparseMatrix"); callGeneric() })            function(x, y = NULL) tcrossprod(x, as(y, "sparseMatrix")))
592    
593    
594  ## FIXME?: In theory, this can be done *FASTER*, in some cases, via tapply1()  ## FIXME?: In theory, this can be done *FASTER*, in some cases, via tapply1()
595  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
596            function(x, y) as(x, "sparseMatrix") %*% y)            function(x, y) as(x, "sparseMatrix") %*% y)
597    setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
598              function(x, y) x %*% as(y, "sparseMatrix"))
599  ## NB: The previous is *not* triggering for  "ddi" o "dgC" (= distance 3)  ## NB: The previous is *not* triggering for  "ddi" o "dgC" (= distance 3)
600  ##     since there's a "ddense" o "Csparse" at dist. 2 => triggers first.  ##     since there's a "ddense" o "Csparse" at dist. 2 => triggers first.
601  ## ==> do this:  ## ==> do this:
# Line 371  Line 607 
607  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*
608  ##       do indeed work by going through sparse (and *not* ddense)!  ##       do indeed work by going through sparse (and *not* ddense)!
609    
 setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),  
           function(x, y) x %*% as(y, "sparseMatrix"))  
610    
611    
612  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),
# Line 395  Line 629 
629  setMethod("solve", signature(a = "diagonalMatrix", b = "Matrix"),  setMethod("solve", signature(a = "diagonalMatrix", b = "Matrix"),
630            solveDiag)            solveDiag)
631    
632    ## Schur()  ---> ./eigen.R
633    
634    
635    
636  ### ---------------- diagonal  o  sparse  -----------------------------  ###---------------- <Ops> (<Arith>, <Logic>, <Compare> ) ----------------------
   
637    
638  ## Use function for several signatures, in order to evade  ## Use function for several signatures, in order to evade
639  ## ambiguous dispatch for "ddi", since there's also Arith(ddense., ddense.)  ## ambiguous dispatch for "ddi", since there's also Arith(ddense., ddense.)
640  diagOdiag <- function(e1,e2) { # result should also be diagonal  diagOdiag <- function(e1,e2) {
641        ## result should also be diagonal _ if possible _
642      r <- callGeneric(.diag.x(e1), .diag.x(e2)) # error if not "compatible"      r <- callGeneric(.diag.x(e1), .diag.x(e2)) # error if not "compatible"
643        ## Check what happens with non-diagonals, i.e. (0 o 0), (FALSE o 0), ...:
644        r00 <- callGeneric(if(is.numeric(e1@x)) 0 else FALSE,
645                           if(is.numeric(e2@x)) 0 else FALSE)
646        if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* diagonal
647      if(is.numeric(r)) {      if(is.numeric(r)) {
648          if(is.numeric(e2@x)) {          if(is.numeric(e2@x)) {
649              e2@x <- r; return(.diag.2N(e2)) }              e2@x <- r; return(.diag.2N(e2)) }
# Line 418  Line 657 
657      e1@x <- r      e1@x <- r
658      .diag.2N(e1)      .diag.2N(e1)
659  }  }
660        else { ## result not diagonal, but at least symmetric:
661            isNum <- (is.numeric(r) || is.numeric(r00))
662            isLog <- (is.logical(r) || is.logical(r00))
663    
664            if(getOption("verbose"))
665                message("exploding  <diag>  o  <diag>  into dense matrix")
666            d <- e1@Dim
667            n <- d[1]
668            stopifnot(length(r) == n)
669            xx <- as.vector(matrix(rbind(r, matrix(r00,n,n)), n,n))
670            newcl <-
671                paste(if(isNum) "d" else if(isLog) {
672                    if(!any(is.na(r)) && !any(is.na(r00))) "n" else "l"
673                } else stop("not yet implemented .. please report")
674                      ,
675                      "syMatrix", sep='')
676    
677            new(newcl, Dim = e1@Dim, Dimnames = e1@Dimnames, x = xx)
678        }
679    }
680    
681    ### This would be *the* way, but we get tons of "ambiguous method dispatch"
682    ## we use this hack instead of signature  x = "diagonalMatrix" :
683    diCls <- names(getClass("diagonalMatrix")@subclasses)
684    if(FALSE) {
685  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "diagonalMatrix"),  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "diagonalMatrix"),
686            diagOdiag)            diagOdiag)
687  ## These two are just for method disambiguation:  } else { ## These are just for method disambiguation:
688  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "diagonalMatrix"),      for(c1 in diCls)
689            diagOdiag)          for(c2 in diCls)
690  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "ddiMatrix"),              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)
691            diagOdiag)  }
692    
693  ## FIXME:    diagonal  o  triangular  |-->  triangular  ## FIXME:    diagonal  o  triangular  |-->  triangular
694  ## -----     diagonal  o  symmetric   |-->  symmetric  ## -----     diagonal  o  symmetric   |-->  symmetric
# Line 435  Line 698 
698  ## For almost everything else, diag* shall be treated "as sparse" :  ## For almost everything else, diag* shall be treated "as sparse" :
699  ## These are cheap implementations via coercion  ## These are cheap implementations via coercion
700    
701  ## for disambiguation  ## For disambiguation --- define this for "sparseMatrix" , then for "ANY";
702  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "sparseMatrix"),  ## and because we can save an .M.kind() call, we use this explicit
703            function(e1,e2) callGeneric(as(e1, "sparseMatrix"), e2))  ## "hack" for all diagonalMatrix *subclasses* instead of just "diagonalMatrix" :
704  setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "diagonalMatrix"),  ##
705            function(e1,e2) callGeneric(e1, as(e2, "sparseMatrix")))  ## ddi*:
706  ## in general:  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "sparseMatrix"),
707  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "ANY"),            function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2))
708            function(e1,e2) callGeneric(as(e1,"sparseMatrix"), e2))  setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "ddiMatrix"),
709  setMethod("Ops", signature(e1 = "ANY", e2 = "diagonalMatrix"),            function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "d")))
710            function(e1,e2) callGeneric(e1, as(e2,"sparseMatrix")))  ## ldi*
711    setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "sparseMatrix"),
712              function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2))
713    setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "ldiMatrix"),
714              function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "l")))
715    
716    ## Ops:  Arith  --> numeric : "dMatrix"
717    ##       Compare --> logical
718    ##       Logic   --> logical: "lMatrix"
719    
720    ##  other = "numeric" : stay diagonal if possible
721    ## ddi*: Arith: result numeric, potentially ddiMatrix
722    setMethod("Arith", signature(e1 = "ddiMatrix", e2 = "numeric"),
723              function(e1,e2) {
724                  n <- e1@Dim[1]; nsq <- n^2
725                  f0 <- callGeneric(0, e2)
726                  if(all(is0(f0))) { # remain diagonal
727                      L1 <- (le <- length(e2)) == 1L
728                      if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)
729                      if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {
730                          e1@diag <- "N"
731                          if(L1) r <- rep.int(r, n)
732                      } else
733                          r <- callGeneric(e1@x, e2)
734                      e1@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
735                      return(e1)
736                  }
737                  callGeneric(diag2tT.u(e1,e2, "d"), e2)
738              })
739    
740    setMethod("Arith", signature(e1 = "numeric", e2 = "ddiMatrix"),
741              function(e1,e2) {
742                  n <- e2@Dim[1]; nsq <- n^2
743                  f0 <- callGeneric(e1, 0)
744                  if(all(is0(f0))) { # remain diagonal
745                      L1 <- (le <- length(e1)) == 1L
746                      if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)
747                      if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {
748                          e2@diag <- "N"
749                          if(L1) r <- rep.int(r, n)
750                      } else
751                          r <- callGeneric(e1, e2@x)
752                      e2@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
753                      return(e2)
754                  }
755                  callGeneric(e1, diag2tT.u(e2,e1, "d"))
756              })
757    
758    ## ldi* Arith --> result numeric, potentially ddiMatrix
759    setMethod("Arith", signature(e1 = "ldiMatrix", e2 = "numeric"),
760              function(e1,e2) {
761                  n <- e1@Dim[1]; nsq <- n^2
762                  f0 <- callGeneric(0, e2)
763                  if(all(is0(f0))) { # remain diagonal
764                      L1 <- (le <- length(e2)) == 1L
765                      if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)
766                      if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {
767                          e1@diag <- "N"
768                          if(L1) r <- rep.int(r, n)
769                      } else
770                          r <- callGeneric(e1@x, e2)
771                      e1 <- copyClass(e1, "ddiMatrix", c("diag", "Dim", "Dimnames"))
772                      e1@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
773                      return(e1)
774                  }
775                  callGeneric(diag2tT.u(e1,e2, "d"), e2)
776              })
777    
778    setMethod("Arith", signature(e1 = "numeric", e2 = "ldiMatrix"),
779              function(e1,e2) {
780                  n <- e2@Dim[1]; nsq <- n^2
781                  f0 <- callGeneric(e1, 0)
782                  if(all(is0(f0))) { # remain diagonal
783                      L1 <- (le <- length(e1)) == 1L
784                      if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)
785                      if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {
786                          e2@diag <- "N"
787                          if(L1) r <- rep.int(r, n)
788                      } else
789                          r <- callGeneric(e1, e2@x)
790                      e2 <- copyClass(e2, "ddiMatrix", c("diag", "Dim", "Dimnames"))
791                      e2@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
792                      return(e2)
793                  }
794                  callGeneric(e1, diag2tT.u(e2,e1, "d"))
795              })
796    
797    ## ddi*: for "Ops" without Arith --> result logical, potentially ldi
798    setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "numeric"),
799              function(e1,e2) {
800                  n <- e1@Dim[1]; nsq <- n^2
801                  f0 <- callGeneric(0, e2)
802                  if(all(is0(f0))) { # remain diagonal
803                      L1 <- (le <- length(e2)) == 1L
804                      if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)
805                      if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {
806                          e1@diag <- "N"
807                          if(L1) r <- rep.int(r, n)
808                      } else
809                          r <- callGeneric(e1@x, e2)
810                      e1 <- copyClass(e1, "ldiMatrix", c("diag", "Dim", "Dimnames"))
811                      e1@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
812                      return(e1)
813                  }
814                  callGeneric(diag2tT.u(e1,e2, "l"), e2)
815              })
816    
817    setMethod("Ops", signature(e1 = "numeric", e2 = "ddiMatrix"),
818              function(e1,e2) {
819                  n <- e2@Dim[1]; nsq <- n^2
820                  f0 <- callGeneric(e1, 0)
821                  if(all(is0(f0))) { # remain diagonal
822                      L1 <- (le <- length(e1)) == 1L
823                      if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)
824                      if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {
825                          e2@diag <- "N"
826                          if(L1) r <- rep.int(r, n)
827                      } else
828                          r <- callGeneric(e1, e2@x)
829                      e2 <- copyClass(e2, "ldiMatrix", c("diag", "Dim", "Dimnames"))
830                      e2@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
831                      return(e2)
832                  }
833                  callGeneric(e1, diag2tT.u(e2,e1, "l"))
834              })
835    
836    ## ldi*: for "Ops" without Arith --> result logical, potentially ldi
837    setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "numeric"),
838              function(e1,e2) {
839                  n <- e1@Dim[1]; nsq <- n^2
840                  f0 <- callGeneric(FALSE, e2)
841                  if(all(is0(f0))) { # remain diagonal
842                      L1 <- (le <- length(e2)) == 1L
843                      if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)
844                      if(e1@diag == "U" && any((r <- callGeneric(TRUE, e2)) != 1)) {
845                          e1@diag <- "N"
846                          if(L1) r <- rep.int(r, n)
847                      } else
848                          r <- callGeneric(e1@x, e2)
849                      e1@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
850                      return(e1)
851                  }
852                  callGeneric(diag2tT.u(e1,e2, "l"), e2)
853              })
854    
855    setMethod("Ops", signature(e1 = "numeric", e2 = "ldiMatrix"),
856              function(e1,e2) {
857                  n <- e2@Dim[1]; nsq <- n^2
858                  f0 <- callGeneric(e1, FALSE)
859                  if(all(is0(f0))) { # remain diagonal
860                      L1 <- (le <- length(e1)) == 1L
861                      if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)
862                      if(e2@diag == "U" && any((r <- callGeneric(e1, TRUE)) != 1)) {
863                          e2@diag <- "N"
864                          if(L1) r <- rep.int(r, n)
865                      } else
866                          r <- callGeneric(e1, e2@x)
867                      e2@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
868                      return(e2)
869                  }
870                  callGeneric(e1, diag2tT.u(e2,e1, "l"))
871              })
872    
873    
874    
875    ## Not {"sparseMatrix", "numeric} :  {"denseMatrix", "matrix", ... }
876    for(other in c("ANY", "Matrix", "dMatrix")) {
877        ## ddi*:
878        setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other),
879                  function(e1,e2) callGeneric(diag2tT.u(e1,e2, "d"), e2))
880        setMethod("Ops", signature(e1 = other, e2 = "ddiMatrix"),
881                  function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "d")))
882        ## ldi*:
883        setMethod("Ops", signature(e1 = "ldiMatrix", e2 = other),
884                  function(e1,e2) callGeneric(diag2tT.u(e1,e2, "l"), e2))
885        setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),
886                  function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "l")))
887    }
888    ## This should *not* dispatch to <dense> methods (in ./Ops.R ), as
889    ##  FALSE & <anything> |-> FALSE : hence result should be diagonal:
890    for(cl in diCls) {
891        setMethod("&", signature(e1 = cl, e2 = "ANY"),
892                  function(e1,e2) e1 & as(e2,"Matrix"))
893        setMethod("&", signature(e1 = "ANY", e2 = cl),
894                  function(e1,e2) as(e1,"Matrix") & e2)
895        for(c2 in c("denseMatrix", "Matrix")) {
896            setMethod("&", signature(e1 = cl, e2 = c2),
897                      function(e1,e2) e1 & Diagonal(x = diag(e2)))
898            setMethod("&", signature(e1 = c2, e2 = cl),
899                      function(e1,e2) Diagonal(x = diag(e1)) & e2)
900        }
901    }
902    
903    
904    ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"
905    ### ----------  any, all: separately here
906    for(cl in diCls) {
907    setMethod("any", cl,
908              function (x, ..., na.rm) {
909                  if(any(x@Dim == 0)) FALSE
910                  else if(x@diag == "U") TRUE else any(x@x, ..., na.rm = na.rm)
911              })
912    setMethod("all",  cl, function (x, ..., na.rm) any(x@Dim == 0))
913    setMethod("prod", cl, function (x, ..., na.rm) as.numeric(any(x@Dim == 0)))
914    
915    setMethod("sum", cl,
916              function(x, ..., na.rm) {
917                  r <- sum(x@x, ..., na.rm = na.rm)# double or integer, correctly
918                  if(x@diag == "U" && !is.na(r)) r + x@Dim[1] else r
919              })
920    }
921    
922    ## The remaining ones are  max, min, range :
923    
924    setMethod("Summary", "ddiMatrix",
925              function(x, ..., na.rm) {
926                  if(any(x@Dim == 0)) callGeneric(numeric(0), ..., na.rm=na.rm)
927                  else if(x@diag == "U")
928                      callGeneric(x@x, 0, 1, ..., na.rm=na.rm)
929                  else callGeneric(x@x, 0, ..., na.rm=na.rm)
930              })
931    setMethod("Summary", "ldiMatrix",
932              function(x, ..., na.rm) {
933                  if(any(x@Dim == 0)) callGeneric(logical(0), ..., na.rm=na.rm)
934                  else if(x@diag == "U")
935                      callGeneric(x@x, FALSE, TRUE, ..., na.rm=na.rm)
936                  else callGeneric(x@x, FALSE, ..., na.rm=na.rm)
937              })
938    
939    
940    
# Line 463  Line 953 
953            function(object) {            function(object) {
954                d <- dim(object)                d <- dim(object)
955                cl <- class(object)                cl <- class(object)
956                cat(sprintf('%d x %d diagonal matrix of class "%s"\n',                cat(sprintf('%d x %d diagonal matrix of class "%s"',
957                            d[1], d[2], cl))                            d[1], d[2], cl))
958                  if(d[1] < 50) {
959                      cat("\n")
960                prDiag(object)                prDiag(object)
961                  } else {
962                      cat(", with diagonal entries\n")
963                      show(diag(object))
964                      invisible(object)
965                  }
966            })            })

Legend:
Removed from v.2052  
changed lines
  Added in v.2209

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