SCM

SCM Repository

[matrix] Diff of /pkg/Matrix/R/diagMatrix.R
ViewVC logotype

Diff of /pkg/Matrix/R/diagMatrix.R

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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

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

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