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

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

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