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

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

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