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 2096, Fri Dec 7 17:44:44 2007 UTC revision 2508, Thu Dec 24 09:47:57 2009 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      ## slightly debatable if we really should return Csparse.. :      r@Dim <- c(i_off[nl+1], j_off[nl + 1])
171      as(ret, "CsparseMatrix")      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    
176  diag2tT <- function(from) {  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    
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      message(sprintf("diagnosing replDiag() -- nargs() = %d\n", nargs()))      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  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "matrix", # 2-col.matrix                       ## 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"),                                  j = "missing", value = "replValue"),
420                   function(x, i, value) {                   function(x,i,j, ..., value) {
421                       if(ncol(i) == 2) {                       if(ncol(i) == 2) {
422                           if(all((ii <- i[,1]) == i[,2])) { # replace in diagonal only                           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                               x@x[ii] <- value
431                               x                               x
432                           } else { ## no longer diagonal, but remain sparse:                           } else { ## no longer diagonal, but remain sparse:
433                               x <- as(x, "sparseMatrix")                               x <- as(x, "TsparseMatrix")
434                               x[i] <- value                               x[i] <- value
435                               x                               x
436                           }                           }
# Line 251  Line 444 
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 262  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")      n <- dimCheck(x,y)[1]
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 305  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")
550      n <- dx[1]      n <- dx[1]
551      as(if(x@diag == "U") y else x@x * y, "Matrix")      as(if(x@diag == "U") y else x@x * y, "Matrix")
552  }  }
   
553  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),
554            diagmatprod)            diagmatprod)
555    ## sneaky .. :
556  formals(diagmatprod) <- alist(x=, y=NULL)  formals(diagmatprod) <- alist(x=, y=NULL)
557  setMethod("crossprod", signature(x = "diagonalMatrix", y = "matrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "matrix"),
558            diagmatprod)            diagmatprod)
 setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "matrix"),  
           diagmatprod)  
559    
560  diagdgeprod <- function(x, y) {  diagGeprod <- function(x, y) {
561      dx <- dim(x)      dx <- dim(x)
562      dy <- dim(y)      dy <- dim(y)
563      if(dx[2] != dy[1]) stop("non-matching dimensions")      if(dx[2] != dy[1]) stop("non-matching dimensions")
# Line 328  Line 565 
565          y@x <- x@x * y@x          y@x <- x@x * y@x
566      y      y
567  }  }
568  setMethod("%*%", signature(x = "diagonalMatrix", y = "dgeMatrix"),  setMethod("%*%", signature(x= "diagonalMatrix", y= "dgeMatrix"), diagGeprod)
569            diagdgeprod, valueClass = "dgeMatrix")  setMethod("%*%", signature(x= "diagonalMatrix", y= "lgeMatrix"), diagGeprod)
570  formals(diagdgeprod) <- alist(x=, y=NULL)  formals(diagGeprod) <- alist(x=, y=NULL)
571  setMethod("crossprod", signature(x = "diagonalMatrix", y = "dgeMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "dgeMatrix"),
572            diagdgeprod, valueClass = "dgeMatrix")            diagGeprod, valueClass = "dgeMatrix")
573    setMethod("crossprod", signature(x = "diagonalMatrix", y = "lgeMatrix"),
574              diagGeprod)
575    
576  setMethod("%*%", signature(x = "matrix", y = "diagonalMatrix"),  matdiagprod <- function(x, y) {
           function(x, y) {  
577                dx <- dim(x)                dx <- dim(x)
578                dy <- dim(y)                dy <- dim(y)
579                if(dx[2] != dy[1]) stop("non-matching dimensions")                if(dx[2] != dy[1]) stop("non-matching dimensions")
580                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]))
581            })  }
582    setMethod("%*%", signature(x = "matrix", y = "diagonalMatrix"),
583              matdiagprod)
584    formals(matdiagprod) <- alist(x=, y=NULL)
585    setMethod("tcrossprod", signature(x = "matrix", y = "diagonalMatrix"),
586              matdiagprod)
587    
588  setMethod("%*%", signature(x = "dgeMatrix", y = "diagonalMatrix"),  gediagprod <- function(x, y) {
           function(x, y) {  
589                dx <- dim(x)                dx <- dim(x)
590                dy <- dim(y)                dy <- dim(y)
591                if(dx[2] != dy[1]) stop("non-matching dimensions")                if(dx[2] != dy[1]) stop("non-matching dimensions")
592                if(y@diag == "N")                if(y@diag == "N")
593                    x@x <- x@x * rep(y@x, each = dx[1])                    x@x <- x@x * rep(y@x, each = dx[1])
594                x                x
595            })  }
596    setMethod("%*%", signature(x= "dgeMatrix", y= "diagonalMatrix"), gediagprod)
597    setMethod("%*%", signature(x= "lgeMatrix", y= "diagonalMatrix"), gediagprod)
598    formals(gediagprod) <- alist(x=, y=NULL)
599    setMethod("tcrossprod", signature(x = "dgeMatrix", y = "diagonalMatrix"),
600              gediagprod)
601    setMethod("tcrossprod", signature(x = "lgeMatrix", y = "diagonalMatrix"),
602              gediagprod)
603    
604  ## crossprod {more of these}  ## crossprod {more of these}
605    
606  ## tcrossprod --- all are not yet there: do the dense ones here:  ## tcrossprod --- all are not yet there: do the dense ones here:
607    
608    setMethod("%*%", signature(x = "diagonalMatrix", y = "denseMatrix"),
609              function(x, y) if(x@diag == "U") y else x %*% as(y, "generalMatrix"))
610    setMethod("%*%", signature(x = "denseMatrix", y = "diagonalMatrix"),
611              function(x, y) if(y@diag == "U") x else as(x, "generalMatrix") %*% y)
612    
613    
614  ## FIXME:  ## FIXME:
615  ## setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "denseMatrix"),  ## setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "denseMatrix"),
616  ##        function(x, y = NULL) {  ##        function(x, y = NULL) {
617  ##           })  ##           })
618    
619  ## setMethod("tcrossprod", signature(x = "denseMatrix", y = "diagonalMatrix"),  Cspdiagprod <- function(x, y) {
620  ##        function(x, y = NULL) {      dx <- dim(x <- .Call(Csparse_diagU2N, x))
621  ##           })      dy <- dim(y)
622        if(dx[2] != dy[1]) stop("non-matching dimensions")
623        ind <- rep.int(seq_len(dx[2]), x@p[-1] - x@p[-dx[2]-1L])
624        if(y@diag == "N")
625            x@x <- x@x * y@x[ind]
626        x
627    }
628    
629    diagCspprod <- function(x, y) {
630        dx <- dim(x)
631        dy <- dim(y <- .Call(Csparse_diagU2N, y))
632        if(dx[2] != dy[1]) stop("non-matching dimensions")
633        if(x@diag == "N")
634            y@x <- y@x * x@x[y@i + 1L]
635        y
636    }
637    
638    setMethod("crossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
639              function(x, y = NULL) diagCspprod(x, y))
640    
641  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
642            function(x, y = NULL) { x <- as(x, "sparseMatrix"); callGeneric() })            function(x, y = NULL) diagCspprod(x, as(y, "CsparseMatrix")))
643    
644    ## Prefer calling diagCspprod to Cspdiagprod if going to transpose anyway
645    ##  x'y == (y'x)'
646    setMethod("crossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
647              function(x, y = NULL) t(diagCspprod(y, x)))
648    
649  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
650            function(x, y = NULL) { y <- as(y, "sparseMatrix"); callGeneric() })            function(x, y = NULL) t(diagCspprod(y, as(x, "Csparsematrix"))))
651    
652    setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
653              function(x, y = NULL) diagCspprod(x, t(y)))
654    
655  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
656            function(x, y = NULL) { x <- as(x, "sparseMatrix"); callGeneric() })            function(x, y = NULL) diagCspprod(x, t(as(y, "CsparseMatrix"))))
657    
658    setMethod("tcrossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
659              function(x, y = NULL) Cspdiagprod(x, y))
660    
661  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
662            function(x, y = NULL) { y <- as(y, "sparseMatrix"); callGeneric() })            function(x, y = NULL) Cspdiagprod(as(x, "CsparseMatrix"), y))
663    
664    setMethod("%*%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
665              function(x, y) diagCspprod(x, y))
666    
 ## FIXME?: In theory, this can be done *FASTER*, in some cases, via tapply1()  
667  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
668            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)!  
669    
670  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
671            function(x, y) x %*% as(y, "sparseMatrix"))            function(x, y) Cspdiagprod(as(x, "CsparseMatrix"), y))
672    
673    setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
674              function(x, y) Cspdiagprod(x, y))
675    
676    ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*
677    ##       do indeed work by going through sparse (and *not* ddense)!
678    
679  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),
680            function(a, b, ...) {            function(a, b, ...) {
# Line 416  Line 696 
696  setMethod("solve", signature(a = "diagonalMatrix", b = "Matrix"),  setMethod("solve", signature(a = "diagonalMatrix", b = "Matrix"),
697            solveDiag)            solveDiag)
698    
699    ## Schur()  ---> ./eigen.R
700    
701    
702    
703  ### ---------------- diagonal  o  sparse  -----------------------------  ###---------------- <Ops> (<Arith>, <Logic>, <Compare> ) ----------------------
   
704    
705  ## Use function for several signatures, in order to evade  ## Use function for several signatures, in order to evade
706  ## ambiguous dispatch for "ddi", since there's also Arith(ddense., ddense.)  ## ambiguous dispatch for "ddi", since there's also Arith(ddense., ddense.)
707  diagOdiag <- function(e1,e2) { # result should also be diagonal  diagOdiag <- function(e1,e2) {
708        ## result should also be diagonal _ if possible _
709      r <- callGeneric(.diag.x(e1), .diag.x(e2)) # error if not "compatible"      r <- callGeneric(.diag.x(e1), .diag.x(e2)) # error if not "compatible"
710        ## Check what happens with non-diagonals, i.e. (0 o 0), (FALSE o 0), ...:
711        r00 <- callGeneric(if(is.numeric(e1@x)) 0 else FALSE,
712                           if(is.numeric(e2@x)) 0 else FALSE)
713        if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* diagonal
714      if(is.numeric(r)) {      if(is.numeric(r)) {
715          if(is.numeric(e2@x)) {          if(is.numeric(e2@x)) {
716              e2@x <- r; return(.diag.2N(e2)) }              e2@x <- r; return(.diag.2N(e2)) }
# Line 439  Line 724 
724      e1@x <- r      e1@x <- r
725      .diag.2N(e1)      .diag.2N(e1)
726  }  }
727        else { ## result not diagonal, but at least symmetric:
728            ## e.g., m == m
729            isNum <- (is.numeric(r) || is.numeric(r00))
730            isLog <- (is.logical(r) || is.logical(r00))
731            Matrix.msg("exploding <diag> o <diag> into dense matrix", .M.level = 2)
732            d <- e1@Dim
733            n <- d[1]
734            stopifnot(length(r) == n)
735            xx <- as.vector(matrix(rbind(r, matrix(r00,n,n)), n,n))
736            newcl <-
737                paste(if(isNum) "d" else if(isLog) {
738                    if(!any(is.na(r)) && !any(is.na(r00))) "n" else "l"
739                } else stop("not yet implemented .. please report")
740                      ,
741                      "syMatrix", sep='')
742    
743            new(newcl, Dim = e1@Dim, Dimnames = e1@Dimnames, x = xx)
744        }
745    }
746    
747    ### This would be *the* way, but we get tons of "ambiguous method dispatch"
748    ## we use this hack instead of signature  x = "diagonalMatrix" :
749    diCls <- names(getClass("diagonalMatrix")@subclasses)
750    if(FALSE) {
751  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "diagonalMatrix"),  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "diagonalMatrix"),
752            diagOdiag)            diagOdiag)
753  ## These two are just for method disambiguation:  } else { ## These are just for method disambiguation:
754  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "diagonalMatrix"),      for(c1 in diCls)
755            diagOdiag)          for(c2 in diCls)
756  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "ddiMatrix"),              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)
757            diagOdiag)  }
758    
759  ## FIXME:    diagonal  o  triangular  |-->  triangular  ## diagonal  o  triangular  |-->  triangular
760  ## -----     diagonal  o  symmetric   |-->  symmetric  ## diagonal  o  symmetric   |-->  symmetric
761  ##    {also when other is sparse: do these "here" --  ##    {also when other is sparse: do these "here" --
762  ##     before conversion to sparse, since that loses "diagonality"}  ##     before conversion to sparse, since that loses "diagonality"}
763    
764  ## For almost everything else, diag* shall be treated "as sparse" :  ## For almost everything else, diag* shall be treated "as sparse" :
765  ## These are cheap implementations via coercion  ## These are cheap implementations via coercion
766    
767  ## for disambiguation  ## For disambiguation --- define this for "sparseMatrix" , then for "ANY";
768  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "sparseMatrix"),  ## and because we can save an .M.kind() call, we use this explicit
769            function(e1,e2) callGeneric(as(e1, "sparseMatrix"), e2))  ## "hack" for all diagonalMatrix *subclasses* instead of just "diagonalMatrix" :
770  setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "diagonalMatrix"),  ##
771            function(e1,e2) callGeneric(e1, as(e2, "sparseMatrix")))  ## ddi*:
772  ## in general:  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "sparseMatrix"),
773  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "ANY"),            function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2))
774            function(e1,e2) callGeneric(as(e1,"sparseMatrix"), e2))  setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "ddiMatrix"),
775  setMethod("Ops", signature(e1 = "ANY", e2 = "diagonalMatrix"),            function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "d")))
776            function(e1,e2) callGeneric(e1, as(e2,"sparseMatrix")))  ## ldi*
777    setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "sparseMatrix"),
778              function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2))
779    setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "ldiMatrix"),
780              function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "l")))
781    
782    ## Ops:  Arith  --> numeric : "dMatrix"
783    ##       Compare --> logical
784    ##       Logic   --> logical: "lMatrix"
785    
786    ##  other = "numeric" : stay diagonal if possible
787    ## ddi*: Arith: result numeric, potentially ddiMatrix
788    setMethod("Arith", signature(e1 = "ddiMatrix", e2 = "numeric"),
789              function(e1,e2) {
790                  n <- e1@Dim[1]; nsq <- n^2
791                  f0 <- callGeneric(0, e2)
792                  if(all(is0(f0))) { # remain diagonal
793                      L1 <- (le <- length(e2)) == 1L
794                      if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)
795                      if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {
796                          e1@diag <- "N"
797                          if(L1) r <- rep.int(r, n)
798                      } else
799                          r <- callGeneric(e1@x, e2)
800                      e1@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
801                      return(e1)
802                  }
803                  callGeneric(diag2tT.u(e1,e2, "d"), e2)
804              })
805    
806    setMethod("Arith", signature(e1 = "numeric", e2 = "ddiMatrix"),
807              function(e1,e2) {
808                  n <- e2@Dim[1]; nsq <- n^2
809                  f0 <- callGeneric(e1, 0)
810                  if(all(is0(f0))) { # remain diagonal
811                      L1 <- (le <- length(e1)) == 1L
812                      if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)
813                      if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {
814                          e2@diag <- "N"
815                          if(L1) r <- rep.int(r, n)
816                      } else
817                          r <- callGeneric(e1, e2@x)
818                      e2@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
819                      return(e2)
820                  }
821                  callGeneric(e1, diag2tT.u(e2,e1, "d"))
822              })
823    
824    ## ldi* Arith --> result numeric, potentially ddiMatrix
825    setMethod("Arith", signature(e1 = "ldiMatrix", e2 = "numeric"),
826              function(e1,e2) {
827                  n <- e1@Dim[1]; nsq <- n^2
828                  f0 <- callGeneric(0, e2)
829                  if(all(is0(f0))) { # remain diagonal
830                      L1 <- (le <- length(e2)) == 1L
831                      if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)
832                      if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {
833                          e1@diag <- "N"
834                          if(L1) r <- rep.int(r, n)
835                      } else
836                          r <- callGeneric(e1@x, e2)
837                      e1 <- copyClass(e1, "ddiMatrix", c("diag", "Dim", "Dimnames"))
838                      e1@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
839                      return(e1)
840                  }
841                  callGeneric(diag2tT.u(e1,e2, "d"), e2)
842              })
843    
844    setMethod("Arith", signature(e1 = "numeric", e2 = "ldiMatrix"),
845              function(e1,e2) {
846                  n <- e2@Dim[1]; nsq <- n^2
847                  f0 <- callGeneric(e1, 0)
848                  if(all(is0(f0))) { # remain diagonal
849                      L1 <- (le <- length(e1)) == 1L
850                      if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)
851                      if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {
852                          e2@diag <- "N"
853                          if(L1) r <- rep.int(r, n)
854                      } else
855                          r <- callGeneric(e1, e2@x)
856                      e2 <- copyClass(e2, "ddiMatrix", c("diag", "Dim", "Dimnames"))
857                      e2@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
858                      return(e2)
859                  }
860                  callGeneric(e1, diag2tT.u(e2,e1, "d"))
861              })
862    
863    ## ddi*: for "Ops" without Arith --> result logical, potentially ldi
864    setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "numeric"),
865              function(e1,e2) {
866                  n <- e1@Dim[1]; nsq <- n^2
867                  f0 <- callGeneric(0, e2)
868                  if(all(is0(f0))) { # remain diagonal
869                      L1 <- (le <- length(e2)) == 1L
870                      if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)
871                      if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {
872                          e1@diag <- "N"
873                          if(L1) r <- rep.int(r, n)
874                      } else
875                          r <- callGeneric(e1@x, e2)
876                      e1 <- copyClass(e1, "ldiMatrix", c("diag", "Dim", "Dimnames"))
877                      e1@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
878                      return(e1)
879                  }
880                  callGeneric(diag2tT.u(e1,e2, "l"), e2)
881              })
882    
883    setMethod("Ops", signature(e1 = "numeric", e2 = "ddiMatrix"),
884              function(e1,e2) {
885                  n <- e2@Dim[1]; nsq <- n^2
886                  f0 <- callGeneric(e1, 0)
887                  if(all(is0(f0))) { # remain diagonal
888                      L1 <- (le <- length(e1)) == 1L
889                      if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)
890                      if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {
891                          e2@diag <- "N"
892                          if(L1) r <- rep.int(r, n)
893                      } else
894                          r <- callGeneric(e1, e2@x)
895                      e2 <- copyClass(e2, "ldiMatrix", c("diag", "Dim", "Dimnames"))
896                      e2@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
897                      return(e2)
898                  }
899                  callGeneric(e1, diag2tT.u(e2,e1, "l"))
900              })
901    
902    ## ldi*: for "Ops" without Arith --> result logical, potentially ldi
903    setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "numeric"),
904              function(e1,e2) {
905                  n <- e1@Dim[1]; nsq <- n^2
906                  f0 <- callGeneric(FALSE, e2)
907                  if(all(is0(f0))) { # remain diagonal
908                      L1 <- (le <- length(e2)) == 1L
909                      if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)
910                      if(e1@diag == "U" && any((r <- callGeneric(TRUE, e2)) != 1)) {
911                          e1@diag <- "N"
912                          if(L1) r <- rep.int(r, n)
913                      } else
914                          r <- callGeneric(e1@x, e2)
915                      e1@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
916                      return(e1)
917                  }
918                  callGeneric(diag2tT.u(e1,e2, "l"), e2)
919              })
920    
921    setMethod("Ops", signature(e1 = "numeric", e2 = "ldiMatrix"),
922              function(e1,e2) {
923                  n <- e2@Dim[1]; nsq <- n^2
924                  f0 <- callGeneric(e1, FALSE)
925                  if(all(is0(f0))) { # remain diagonal
926                      L1 <- (le <- length(e1)) == 1L
927                      if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)
928                      if(e2@diag == "U" && any((r <- callGeneric(e1, TRUE)) != 1)) {
929                          e2@diag <- "N"
930                          if(L1) r <- rep.int(r, n)
931                      } else
932                          r <- callGeneric(e1, e2@x)
933                      e2@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
934                      return(e2)
935                  }
936                  callGeneric(e1, diag2tT.u(e2,e1, "l"))
937              })
938    
939    
940    
941    ## Not {"sparseMatrix", "numeric} :  {"denseMatrix", "matrix", ... }
942    for(other in c("ANY", "Matrix", "dMatrix")) {
943        ## ddi*:
944        setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other),
945                  function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2))
946        setMethod("Ops", signature(e1 = other, e2 = "ddiMatrix"),
947                  function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "d")))
948        ## ldi*:
949        setMethod("Ops", signature(e1 = "ldiMatrix", e2 = other),
950                  function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2))
951        setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),
952                  function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "l")))
953    }
954    
955    ## Direct subclasses of "denseMatrix": currently ddenseMatrix, ldense... :
956    dense.subCl <- local({ dM.scl <- getClass("denseMatrix")@subclasses
957                           names(dM.scl)[sapply(dM.scl, slot, "distance") == 1] })
958    for(DI in diCls) {
959        for(c2 in c(dense.subCl, "Matrix")) {
960            for(Fun in c("*", "^", "&")) {
961                setMethod(Fun, signature(e1 = DI, e2 = c2),
962                          function(e1,e2) callGeneric(e1, Diagonal(x = diag(e2))))
963                setMethod(Fun, signature(e1 = c2, e2 = DI),
964                          function(e1,e2) callGeneric(Diagonal(x = diag(e1)), e2))
965            }
966            ## NB: This arguably implicitly uses  0/0 :== 0  to keep diagonality
967            for(Fun in c("%%", "%/%", "/")) {
968                setMethod(Fun, signature(e1 = DI, e2 = c2),
969                          function(e1,e2) callGeneric(e1, Diagonal(x = diag(e2))))
970            }
971        }
972    }
973    
974    
975    ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"
976    ### ----------   the last 4: separately here
977    for(cl in diCls) {
978    setMethod("any", cl,
979              function (x, ..., na.rm) {
980                  if(any(x@Dim == 0)) FALSE
981                  else if(x@diag == "U") TRUE else any(x@x, ..., na.rm = na.rm)
982              })
983    setMethod("all",  cl, function (x, ..., na.rm) {
984        n <- x@Dim[1]
985        if(n >= 2) FALSE
986        else if(n == 0 || x@diag == "U") TRUE
987        else all(x@x, ..., na.rm = na.rm)
988    })
989    setMethod("prod", cl, function (x, ..., na.rm) {
990        n <- x@Dim[1]
991        if(n >= 2) 0
992        else if(n == 0 || x@diag == "U") 1
993        else ## n == 1, diag = "N" :
994            prod(x@x, ..., na.rm = na.rm)
995    })
996    
997    setMethod("sum", cl,
998              function(x, ..., na.rm) {
999                  r <- sum(x@x, ..., na.rm = na.rm)# double or integer, correctly
1000                  if(x@diag == "U" && !is.na(r)) r + x@Dim[1] else r
1001              })
1002    }
1003    
1004    ## The remaining ones are  max, min, range :
1005    
1006    setMethod("Summary", "ddiMatrix",
1007              function(x, ..., na.rm) {
1008                  if(any(x@Dim == 0)) callGeneric(numeric(0), ..., na.rm=na.rm)
1009                  else if(x@diag == "U")
1010                      callGeneric(x@x, 0, 1, ..., na.rm=na.rm)
1011                  else callGeneric(x@x, 0, ..., na.rm=na.rm)
1012              })
1013    setMethod("Summary", "ldiMatrix",
1014              function(x, ..., na.rm) {
1015                  if(any(x@Dim == 0)) callGeneric(logical(0), ..., na.rm=na.rm)
1016                  else if(x@diag == "U")
1017                      callGeneric(x@x, FALSE, TRUE, ..., na.rm=na.rm)
1018                  else callGeneric(x@x, FALSE, ..., na.rm=na.rm)
1019              })
1020    
1021    
1022    
# Line 480  Line 1031 
1031      invisible(x)      invisible(x)
1032  }  }
1033    
1034    ## somewhat consistent with "print" for sparseMatrix :
1035    setMethod("print", signature(x = "diagonalMatrix"), prDiag)
1036    
1037  setMethod("show", signature(object = "diagonalMatrix"),  setMethod("show", signature(object = "diagonalMatrix"),
1038            function(object) {            function(object) {
1039                d <- dim(object)                d <- dim(object)
1040                cl <- class(object)                cl <- class(object)
1041                cat(sprintf('%d x %d diagonal matrix of class "%s"\n',                cat(sprintf('%d x %d diagonal matrix of class "%s"',
1042                            d[1], d[2], cl))                            d[1], d[2], cl))
1043                  if(d[1] < 50) {
1044                      cat("\n")
1045                prDiag(object)                prDiag(object)
1046                  } else {
1047                      cat(", with diagonal entries\n")
1048                      show(diag(object))
1049                      invisible(object)
1050                  }
1051            })            })

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

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