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 2157, Tue Mar 25 15:00:01 2008 UTC pkg/Matrix/R/diagMatrix.R revision 2608, Tue Aug 10 08:23:48 2010 UTC
# Line 32  Line 32 
32      }      }
33  }  }
34    
35  ## Pkg 'spdep' had (relatively slow) versions of this as_dsCMatrix_I()  .sparseDiagonal <- function(n, x = rep.int(1,m), uplo = "U",
36  .symDiagonal <- function(n, x = rep.int(1,n), uplo = "U") {                              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)      stopifnot(n == (n. <- as.integer(n)), (n <- n.) >= 0)
40      if((lx <- length(x)) == 1) x <- rep.int(x, n)      if(!missing(cols))
41      else if(lx != n) stop("length(x) must be 1 or n")          stopifnot(0 <= (cols <- as.integer(cols)), cols < n)
42      cls <-      m <- length(cols)
43          if(is.double(x)) "dsCMatrix"      if(missing(kind))
44          else if(is.logical(x)) "lsCMatrix"          kind <-
45                if(is.double(x)) "d"
46                else if(is.logical(x)) "l"
47          else { ## for now          else { ## for now
48              storage.mode(x) <- "double"              storage.mode(x) <- "double"
49              "dsCMatrix"                  "d"
50          }          }
51      new(cls, Dim = c(n,n), x = x, uplo = uplo,      else stopifnot(any(kind == c("d","l","n")))
52          i = if(n) 0:(n - 1L) else integer(0), p = 0:n)      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      ## slightly debatable if we really should return Csparse.. :      else if(all(tri) && { ULs <- sapply(Tlst, slot, "uplo")##  "U" or "L"
149      as(ret, "CsparseMatrix")                            all(ULs[1L] == ULs[-1L]) } ## all upper or all lower
150           ){ ## *triangular* result
151    
152            if(use.n) { ## return nsparseMatrix :
153                r <- new("ntTMatrix")
154            } else {
155                r <- new("dtTMatrix")
156                r@x <- unlist(lapply(Tlst, slot, "x"))
157            }
158            r@uplo <- ULs[1L]
159        }
160        else {
161            if(any(sym))
162                Tlst[sym] <- lapply(Tlst[sym], as, "generalMatrix")
163            if(use.n) { ## return nsparseMatrix :
164                r <- new("ngTMatrix")
165            } else {
166                r <- new("dgTMatrix")
167                r@x <- unlist(lapply(Tlst, slot, "x"))
168            }
169        }
170        r@Dim <- c(i_off[nl+1], j_off[nl + 1])
171        r@i <- unlist(lapply(1:nl, function(k) Tlst[[k]]@i + i_off[k]))
172        r@j <- unlist(lapply(1:nl, function(k) Tlst[[k]]@j + j_off[k]))
173        r
174    }
175    
176    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    
# Line 107  Line 216 
216  diag2tT.u <- function(d, x, kind = .M.kind(d))  diag2tT.u <- function(d, x, kind = .M.kind(d))
217      .diag2tT(d, uplo = if(is(x,"triangularMatrix")) x@uplo else "U", kind)      .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    
229  ## In order to evade method dispatch ambiguity warnings,  ## In order to evade method dispatch ambiguity warnings,
230  ## and because we can save a .M.kind() call, we use this explicit  ## and because we can save a .M.kind() call, we use this explicit
# Line 115  Line 233 
233  ## ddi*:  ## ddi*:
234  diag2tT <- function(from) .diag2tT(from, "U", "d")  diag2tT <- function(from) .diag2tT(from, "U", "d")
235  setAs("ddiMatrix", "triangularMatrix", diag2tT)  setAs("ddiMatrix", "triangularMatrix", diag2tT)
236  setAs("ddiMatrix", "sparseMatrix", diag2tT)  ##_no_longer_ setAs("ddiMatrix", "sparseMatrix", diag2tT)
237  ## needed too (otherwise <dense> -> Tsparse is taken):  ## needed too (otherwise <dense> -> Tsparse is taken):
238  setAs("ddiMatrix", "TsparseMatrix", diag2tT)  setAs("ddiMatrix", "TsparseMatrix", diag2tT)
239    setAs("ddiMatrix", "dsparseMatrix", diag2tT)
240  setAs("ddiMatrix", "CsparseMatrix",  setAs("ddiMatrix", "CsparseMatrix",
241        function(from) as(.diag2tT(from, "U", "d"), "CsparseMatrix"))        function(from) as(.diag2tT(from, "U", "d"), "CsparseMatrix"))
242  setAs("ddiMatrix", "symmetricMatrix",  setAs("ddiMatrix", "symmetricMatrix",
# Line 126  Line 245 
245  ## ldi*:  ## ldi*:
246  diag2tT <- function(from) .diag2tT(from, "U", "l")  diag2tT <- function(from) .diag2tT(from, "U", "l")
247  setAs("ldiMatrix", "triangularMatrix", diag2tT)  setAs("ldiMatrix", "triangularMatrix", diag2tT)
248  setAs("ldiMatrix", "sparseMatrix", diag2tT)  ##_no_longer_ setAs("ldiMatrix", "sparseMatrix", diag2tT)
249  ## needed too (otherwise <dense> -> Tsparse is taken):  ## needed too (otherwise <dense> -> Tsparse is taken):
250  setAs("ldiMatrix", "TsparseMatrix", diag2tT)  setAs("ldiMatrix", "TsparseMatrix", diag2tT)
251    setAs("ldiMatrix", "lsparseMatrix", diag2tT)
252  setAs("ldiMatrix", "CsparseMatrix",  setAs("ldiMatrix", "CsparseMatrix",
253        function(from) as(.diag2tT(from, "U", "l"), "CsparseMatrix"))        function(from) as(.diag2tT(from, "U", "l"), "CsparseMatrix"))
254  setAs("ldiMatrix", "symmetricMatrix",  setAs("ldiMatrix", "symmetricMatrix",
# Line 143  Line 263 
263                Dim = from@Dim, Dimnames = from@Dimnames)                Dim = from@Dim, Dimnames = from@Dimnames)
264        })        })
265    
266    setAs("diagonalMatrix", "nsparseMatrix", function(from) as(from, "nMatrix"))
267    
268    ## Cheap fast substitute for diag() which *does* preserve the mode of x :
269    mkDiag <- function(x, n) {
270        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  setAs("diagonalMatrix", "matrix",  setAs("diagonalMatrix", "matrix",
276        function(from) {        function(from) {
277            n <- from@Dim[1]            ## want "ldiMatrix" -> <logical> "matrix" :
278            diag(x = if(from@diag == "U") { if(is.numeric(from@x)) 1. else TRUE            mkDiag(if(from@diag == "U") as1(from@x) else from@x,
279                                       } else from@x,                   n = from@Dim[1])
                nrow = n, ncol = n)  
280        })        })
281    
282  setMethod("as.vector", signature(x = "diagonalMatrix", mode="missing"),  setMethod("as.vector", signature(x = "diagonalMatrix", mode="missing"),
283            function(x, mode) {            function(x, mode) {
284                n <- x@Dim[1]                n <- x@Dim[1]
285                mod <- mode(x@x)                mod.x <- mode(x@x)
286                r <- vector(mod, length = n^2)                r <- vector(mod.x, length = n^2)
287                if(n)                if(n)
288                    r[1 + 0:(n - 1) * (n + 1)] <-                    r[1 + 0:(n - 1) * (n + 1)] <-
289                        if(x@diag == "U")                        if(x@diag == "U") as1(mod=mod.x) else x@x
                           switch(mod, "integer"= 1L,  
                                  "numeric"= 1, "logical"= TRUE)  
                       else x@x  
290                r                r
291            })            })
292    
293  setAs("diagonalMatrix", "generalMatrix", # prefer sparse:  setAs("diagonalMatrix", "generalMatrix", # prefer sparse:
294        function(from) as(as(from, "CsparseMatrix"), "generalMatrix"))        function(from) as(as(from, "CsparseMatrix"), "generalMatrix"))
295    
296  .diag.x <- function(m) {  setAs("diagonalMatrix", "denseMatrix",
297      if(m@diag == "U")        function(from) as(as(from, "CsparseMatrix"), "denseMatrix"))
298          rep.int(if(is.numeric(m@x)) 1. else TRUE, m@Dim[1])  
299      else m@x  .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) {  .diag.2N <- function(m) {
302      if(m@diag == "U") m@diag <- "N"      if(m@diag == "U") m@diag <- "N"
303      m      m
304  }  }
305    
 if(FALSE) {  
 ## given the above, the following  4  coercions should be all unneeded;  
 ## we prefer triangular to general:  
 setAs("ddiMatrix", "dgTMatrix",  
       function(from) {  
           .Deprecated("as(, \"sparseMatrix\")")  
           n <- from@Dim[1]  
           i <- seq_len(n) - 1L  
           new("dgTMatrix", i = i, j = i, x = .diag.x(from),  
               Dim = c(n,n), Dimnames = from@Dimnames) })  
   
 setAs("ddiMatrix", "dgCMatrix",  
       function(from) as(as(from, "sparseMatrix"), "dgCMatrix"))  
   
 setAs("ldiMatrix", "lgTMatrix",  
       function(from) {  
           .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) })  
   
 setAs("ldiMatrix", "lgCMatrix",  
       function(from) as(as(from, "lgTMatrix"), "lgCMatrix"))  
 }  
   
   
 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 247  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        })        })
352    
353    
354  ## In order to evade method dispatch ambiguity warnings,  setMethod("diag", signature(x = "diagonalMatrix"),
 ## we use this hack instead of signature  x = "diagonalMatrix" :  
 diCls <- names(getClass("diagonalMatrix")@subclasses)  
 for(cls in diCls) {  
     setMethod("diag", signature(x = cls),  
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) {  subDiag <- function(x, i, j, ..., drop) {
358      x <- as(x, "sparseMatrix")      x <- as(x, "CsparseMatrix") ## << was "TsparseMatrix" (Csparse is faster now)
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(isS4(x) && isDiagonal(x)) as(x, "diagonalMatrix") else x      if(isS4(x) && isDiagonal(x)) as(x, "diagonalMatrix") else x
# Line 282  Line 369 
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, j, ..., 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, i, j, ..., drop) subDiag(x, j=j, drop=drop))            function(x, i, j, ..., drop) subDiag(x, j=j, drop=drop))
# Line 292  Line 385 
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  ## Only(?) current bug:  x[i] <- value  is wrong when  i is *vector*  ## Only(?) current bug:  x[i] <- value  is wrong when  i is *vector*
387  replDiag <- function(x, i, j, ..., value) {  replDiag <- function(x, i, j, ..., value) {
388      x <- as(x, "sparseMatrix")      x <- as(x, "TsparseMatrix")
389      if(missing(i))      if(missing(i))
390          x[, j] <- value          x[, j] <- value
391      else if(missing(j)) { ##  x[i , ] <- v  *OR*   x[i] <- v      else if(missing(j)) { ##  x[i , ] <- v  *OR*   x[i] <- v
# Line 321  Line 414 
414                           replDiag(x, i=i, , value=value)                           replDiag(x, i=i, , value=value)
415                   })                   })
416    
417  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "matrix", # 2-col.matrix  setReplaceMethod("[", signature(x = "diagonalMatrix",
418                                    i = "matrix", # 2-col.matrix
419                                  j = "missing", value = "replValue"),                                  j = "missing", value = "replValue"),
420                   function(x,i,j, ..., 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 345  Line 446 
446                                  j = "index", value = "replValue"),                                  j = "index", value = "replValue"),
447                   function(x,i,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"),
474            function(x) { x@Dimnames <- x@Dimnames[2:1] ; x })            function(x) { x@Dimnames <- x@Dimnames[2:1] ; x })
# Line 361  Line 485 
485    
486  setMethod("chol", signature(x = "ddiMatrix"),  setMethod("chol", signature(x = "ddiMatrix"),
487            function(x, pivot, ...) {            function(x, pivot, ...) {
488                  if(x@diag == "U") return(x)
489                  ## else
490                if(any(x@x < 0))                if(any(x@x < 0))
491                    stop("chol() is undefined for diagonal matrix with negative entries")                    stop("chol() is undefined for diagonal matrix with negative entries")
492                x@x <- sqrt(x@x)                x@x <- sqrt(x@x)
# Line 369  Line 495 
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 401  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 424  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        if(is(x, "compMatrix") && length(xf <- x@factors)) {
627            ## instead of dropping all factors, be smart about some
628            ## TODO ......
629            x@factors <- list()
630        }
631        x
632    }
633    
634    diagCspprod <- function(x, y) {
635        dx <- dim(x)
636        dy <- dim(y <- .Call(Csparse_diagU2N, y))
637        if(dx[2] != dy[1]) stop("non-matching dimensions")
638        if(x@diag == "N")
639            y@x <- y@x * x@x[y@i + 1L]
640        if(is(y, "compMatrix") && length(yf <- y@factors)) {
641            ## instead of dropping all factors, be smart about some
642            ## TODO
643            keep <- character()
644            if(iLU <- names(yf) == "LU") {
645                ## TODO keep <- "LU"
646            }
647            y@factors <- yf[keep]
648        }
649        y
650    }
651    
652    setMethod("crossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
653              function(x, y = NULL) diagCspprod(x, y))
654    
655  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
656            function(x, y = NULL) crossprod(as(x, "sparseMatrix"), y))            function(x, y = NULL) diagCspprod(x, as(y, "CsparseMatrix")))
657    
658    ## Prefer calling diagCspprod to Cspdiagprod if going to transpose anyway
659    ##  x'y == (y'x)'
660    setMethod("crossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
661              function(x, y = NULL) t(diagCspprod(y, x)))
662    
663  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
664            function(x, y = NULL) crossprod(x, as(y, "sparseMatrix")))            function(x, y = NULL) t(diagCspprod(y, as(x, "Csparsematrix"))))
665    
666    setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
667              function(x, y = NULL) diagCspprod(x, t(y)))
668    
669  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
670            function(x, y = NULL) tcrossprod(as(x, "sparseMatrix"), y))            function(x, y = NULL) diagCspprod(x, t(as(y, "CsparseMatrix"))))
671    
672    setMethod("tcrossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
673              function(x, y = NULL) Cspdiagprod(x, y))
674    
675  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
676            function(x, y = NULL) tcrossprod(x, as(y, "sparseMatrix")))            function(x, y = NULL) Cspdiagprod(as(x, "CsparseMatrix"), y))
677    
678    setMethod("%*%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
679              function(x, y) diagCspprod(x, y))
680    
 ## FIXME?: In theory, this can be done *FASTER*, in some cases, via tapply1()  
681  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
682            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)!  
683    
684  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
685            function(x, y) x %*% as(y, "sparseMatrix"))            function(x, y) Cspdiagprod(as(x, "CsparseMatrix"), y))
686    
687    setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
688              function(x, y) Cspdiagprod(x, y))
689    
690    ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*
691    ##       do indeed work by going through sparse (and *not* ddense)!
692    
693  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),
694            function(a, b, ...) {            function(a, b, ...) {
# Line 516  Line 714 
714    
715    
716    
717  ### ---------------- diagonal  o  sparse  -----------------------------  ###---------------- <Ops> (<Arith>, <Logic>, <Compare> ) ----------------------
   
718    
719  ## Use function for several signatures, in order to evade  ## Use function for several signatures, in order to evade
720  ## ambiguous dispatch for "ddi", since there's also Arith(ddense., ddense.)  ## ambiguous dispatch for "ddi", since there's also Arith(ddense., ddense.)
721  diagOdiag <- function(e1,e2) { # result should also be diagonal  diagOdiag <- function(e1,e2) {
722        ## result should also be diagonal _ if possible _
723      r <- callGeneric(.diag.x(e1), .diag.x(e2)) # error if not "compatible"      r <- callGeneric(.diag.x(e1), .diag.x(e2)) # error if not "compatible"
724        ## Check what happens with non-diagonals, i.e. (0 o 0), (FALSE o 0), ...:
725        r00 <- callGeneric(if(is.numeric(e1@x)) 0 else FALSE,
726                           if(is.numeric(e2@x)) 0 else FALSE)
727        if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* diagonal
728      if(is.numeric(r)) {      if(is.numeric(r)) {
729          if(is.numeric(e2@x)) {          if(is.numeric(e2@x)) {
730              e2@x <- r; return(.diag.2N(e2)) }              e2@x <- r; return(.diag.2N(e2)) }
# Line 536  Line 738 
738      e1@x <- r      e1@x <- r
739      .diag.2N(e1)      .diag.2N(e1)
740  }  }
741        else { ## result not diagonal, but at least symmetric:
742            ## e.g., m == m
743            isNum <- (is.numeric(r) || is.numeric(r00))
744            isLog <- (is.logical(r) || is.logical(r00))
745            Matrix.msg("exploding <diag> o <diag> into dense matrix", .M.level = 2)
746            d <- e1@Dim
747            n <- d[1]
748            stopifnot(length(r) == n)
749            xx <- as.vector(matrix(rbind(r, matrix(r00,n,n)), n,n))
750            newcl <-
751                paste(if(isNum) "d" else if(isLog) {
752                    if(!any(is.na(r)) && !any(is.na(r00))) "n" else "l"
753                } else stop("not yet implemented .. please report")
754                      ,
755                      "syMatrix", sep='')
756    
757            new(newcl, Dim = e1@Dim, Dimnames = e1@Dimnames, x = xx)
758        }
759    }
760    
761    ### This would be *the* way, but we get tons of "ambiguous method dispatch"
762    ## we use this hack instead of signature  x = "diagonalMatrix" :
763    diCls <- names(getClass("diagonalMatrix")@subclasses)
764    if(FALSE) {
765  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "diagonalMatrix"),  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "diagonalMatrix"),
766            diagOdiag)            diagOdiag)
767  ## These two are just for method disambiguation:  } else { ## These are just for method disambiguation:
768  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "diagonalMatrix"),      for(c1 in diCls)
769            diagOdiag)          for(c2 in diCls)
770  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "ddiMatrix"),              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)
771            diagOdiag)  }
772    
773  ## FIXME:    diagonal  o  triangular  |-->  triangular  ## diagonal  o  triangular  |-->  triangular
774  ## -----     diagonal  o  symmetric   |-->  symmetric  ## diagonal  o  symmetric   |-->  symmetric
775  ##    {also when other is sparse: do these "here" --  ##    {also when other is sparse: do these "here" --
776  ##     before conversion to sparse, since that loses "diagonality"}  ##     before conversion to sparse, since that loses "diagonality"}
777    
# Line 559  Line 784 
784  ##  ##
785  ## ddi*:  ## ddi*:
786  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "sparseMatrix"),  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "sparseMatrix"),
787            function(e1,e2) callGeneric(diag2tT.u(e1,e2, "d"), e2))            function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2))
788  setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "ddiMatrix"),  setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "ddiMatrix"),
789            function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "d")))            function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "d")))
790  ## ldi*  ## ldi*
791  setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "sparseMatrix"),  setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "sparseMatrix"),
792            function(e1,e2) callGeneric(diag2tT.u(e1,e2, "l"), e2))            function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2))
793  setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "ldiMatrix"),  setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "ldiMatrix"),
794            function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "l")))            function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "l")))
795    
796    ## Ops:  Arith  --> numeric : "dMatrix"
797    ##       Compare --> logical
798    ##       Logic   --> logical: "lMatrix"
799    
800  ##  other = "numeric" : stay diagonal if possible  ##  other = "numeric" : stay diagonal if possible
801  ## ddi*:  ## ddi*: Arith: result numeric, potentially ddiMatrix
802  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "numeric"),  setMethod("Arith", signature(e1 = "ddiMatrix", e2 = "numeric"),
803            function(e1,e2) {            function(e1,e2) {
804                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]; nsq <- n^2
805                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
806                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
807                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
# Line 582  Line 811 
811                        if(L1) r <- rep.int(r, n)                        if(L1) r <- rep.int(r, n)
812                    } else                    } else
813                        r <- callGeneric(e1@x, e2)                        r <- callGeneric(e1@x, e2)
814                    e1@x <- if(L1) r else r[1L + n*(0:(n-1L))]                    e1@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
815                    return(e1)                    return(e1)
816                }                }
817                callGeneric(diag2tT.u(e1,e2, "d"), e2)                callGeneric(diag2tT.u(e1,e2, "d"), e2)
818            })            })
819    
820  setMethod("Ops", signature(e1 = "numeric", e2 = "ddiMatrix"),  setMethod("Arith", signature(e1 = "numeric", e2 = "ddiMatrix"),
821            function(e1,e2) {            function(e1,e2) {
822                n <- e2@Dim[1]; nsq <- n*n                n <- e2@Dim[1]; nsq <- n^2
823                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
824                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
825                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
# Line 600  Line 829 
829                        if(L1) r <- rep.int(r, n)                        if(L1) r <- rep.int(r, n)
830                    } else                    } else
831                        r <- callGeneric(e1, e2@x)                        r <- callGeneric(e1, e2@x)
832                    e2@x <- if(L1) r else r[1L + n*(0:(n-1L))]                    e2@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
833                    return(e2)                    return(e2)
834                }                }
835                callGeneric(e1, diag2tT.u(e2,e1, "d"))                callGeneric(e1, diag2tT.u(e2,e1, "d"))
836            })            })
837  ## ldi*:  
838    ## ldi* Arith --> result numeric, potentially ddiMatrix
839    setMethod("Arith", signature(e1 = "ldiMatrix", e2 = "numeric"),
840              function(e1,e2) {
841                  n <- e1@Dim[1]; nsq <- n^2
842                  f0 <- callGeneric(0, e2)
843                  if(all(is0(f0))) { # remain diagonal
844                      L1 <- (le <- length(e2)) == 1L
845                      if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)
846                      if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {
847                          e1@diag <- "N"
848                          if(L1) r <- rep.int(r, n)
849                      } else
850                          r <- callGeneric(e1@x, e2)
851                      e1 <- copyClass(e1, "ddiMatrix", c("diag", "Dim", "Dimnames"))
852                      e1@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
853                      return(e1)
854                  }
855                  callGeneric(diag2tT.u(e1,e2, "d"), e2)
856              })
857    
858    setMethod("Arith", signature(e1 = "numeric", e2 = "ldiMatrix"),
859              function(e1,e2) {
860                  n <- e2@Dim[1]; nsq <- n^2
861                  f0 <- callGeneric(e1, 0)
862                  if(all(is0(f0))) { # remain diagonal
863                      L1 <- (le <- length(e1)) == 1L
864                      if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)
865                      if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {
866                          e2@diag <- "N"
867                          if(L1) r <- rep.int(r, n)
868                      } else
869                          r <- callGeneric(e1, e2@x)
870                      e2 <- copyClass(e2, "ddiMatrix", c("diag", "Dim", "Dimnames"))
871                      e2@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
872                      return(e2)
873                  }
874                  callGeneric(e1, diag2tT.u(e2,e1, "d"))
875              })
876    
877    ## ddi*: for "Ops" without Arith --> result logical, potentially ldi
878    setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "numeric"),
879              function(e1,e2) {
880                  n <- e1@Dim[1]; nsq <- n^2
881                  f0 <- callGeneric(0, e2)
882                  if(all(is0(f0))) { # remain diagonal
883                      L1 <- (le <- length(e2)) == 1L
884                      if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)
885                      if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {
886                          e1@diag <- "N"
887                          if(L1) r <- rep.int(r, n)
888                      } else
889                          r <- callGeneric(e1@x, e2)
890                      e1 <- copyClass(e1, "ldiMatrix", c("diag", "Dim", "Dimnames"))
891                      e1@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
892                      return(e1)
893                  }
894                  callGeneric(diag2tT.u(e1,e2, "l"), e2)
895              })
896    
897    setMethod("Ops", signature(e1 = "numeric", e2 = "ddiMatrix"),
898              function(e1,e2) {
899                  n <- e2@Dim[1]; nsq <- n^2
900                  f0 <- callGeneric(e1, 0)
901                  if(all(is0(f0))) { # remain diagonal
902                      L1 <- (le <- length(e1)) == 1L
903                      if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)
904                      if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {
905                          e2@diag <- "N"
906                          if(L1) r <- rep.int(r, n)
907                      } else
908                          r <- callGeneric(e1, e2@x)
909                      e2 <- copyClass(e2, "ldiMatrix", c("diag", "Dim", "Dimnames"))
910                      e2@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
911                      return(e2)
912                  }
913                  callGeneric(e1, diag2tT.u(e2,e1, "l"))
914              })
915    
916    ## ldi*: for "Ops" without Arith --> result logical, potentially ldi
917  setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "numeric"),  setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "numeric"),
918            function(e1,e2) {            function(e1,e2) {
919                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]; nsq <- n^2
920                f0 <- callGeneric(FALSE, e2)                f0 <- callGeneric(FALSE, e2)
921                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
922                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
# Line 618  Line 926 
926                        if(L1) r <- rep.int(r, n)                        if(L1) r <- rep.int(r, n)
927                    } else                    } else
928                        r <- callGeneric(e1@x, e2)                        r <- callGeneric(e1@x, e2)
929                    e1@x <- if(L1) r else r[1L + n*(0:(n-1L))]                    e1@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
930                    return(e1)                    return(e1)
931                }                }
932                callGeneric(diag2tT.u(e1,e2, "l"), e2)                callGeneric(diag2tT.u(e1,e2, "l"), e2)
# Line 626  Line 934 
934    
935  setMethod("Ops", signature(e1 = "numeric", e2 = "ldiMatrix"),  setMethod("Ops", signature(e1 = "numeric", e2 = "ldiMatrix"),
936            function(e1,e2) {            function(e1,e2) {
937                n <- e2@Dim[1]; nsq <- n*n                n <- e2@Dim[1]; nsq <- n^2
938                f0 <- callGeneric(e1, FALSE)                f0 <- callGeneric(e1, FALSE)
939                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
940                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
# Line 636  Line 944 
944                        if(L1) r <- rep.int(r, n)                        if(L1) r <- rep.int(r, n)
945                    } else                    } else
946                        r <- callGeneric(e1, e2@x)                        r <- callGeneric(e1, e2@x)
947                    e2@x <- if(L1) r else r[1L + n*(0:(n-1L))]                    e2@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
948                    return(e2)                    return(e2)
949                }                }
950                callGeneric(e1, diag2tT.u(e2,e1, "l"))                callGeneric(e1, diag2tT.u(e2,e1, "l"))
# Line 645  Line 953 
953    
954    
955  ## Not {"sparseMatrix", "numeric} :  {"denseMatrix", "matrix", ... }  ## Not {"sparseMatrix", "numeric} :  {"denseMatrix", "matrix", ... }
956    for(other in c("ANY", "Matrix", "dMatrix")) {
957  ## ddi*:  ## ddi*:
958  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "ANY"),      setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other),
959            function(e1,e2) callGeneric(diag2tT.u(e1,e2, "d"), e2))                function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2))
960  setMethod("Ops", signature(e1 = "ANY", e2 = "ddiMatrix"),      setMethod("Ops", signature(e1 = other, e2 = "ddiMatrix"),
961            function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "d")))                function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "d")))
962  ## ldi*:  ## ldi*:
963  setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "ANY"),      setMethod("Ops", signature(e1 = "ldiMatrix", e2 = other),
964            function(e1,e2) callGeneric(diag2tT.u(e1,e2, "l"), e2))                function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2))
965  setMethod("Ops", signature(e1 = "ANY", e2 = "ldiMatrix"),      setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),
966            function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "l")))                function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "l")))
967    }
968    
969    ## Direct subclasses of "denseMatrix": currently ddenseMatrix, ldense... :
970    if(FALSE)## too general, would contain  denseModelMatrix:
971    dense.subCl <- local({ dM.scl <- getClass("denseMatrix")@subclasses
972                           names(dM.scl)[sapply(dM.scl, slot, "distance") == 1] })
973    dense.subCl <- paste(c("d","l","n"), "denseMatrix", sep="")
974    for(DI in diCls) {
975        for(c2 in c(dense.subCl, "Matrix")) {
976            for(Fun in c("*", "^", "&")) {
977                setMethod(Fun, signature(e1 = DI, e2 = c2),
978                          function(e1,e2) callGeneric(e1, Diagonal(x = diag(e2))))
979                setMethod(Fun, signature(e1 = c2, e2 = DI),
980                          function(e1,e2) callGeneric(Diagonal(x = diag(e1)), e2))
981            }
982            ## NB: This arguably implicitly uses  0/0 :== 0  to keep diagonality
983            for(Fun in c("%%", "%/%", "/")) {
984                setMethod(Fun, signature(e1 = DI, e2 = c2),
985                          function(e1,e2) callGeneric(e1, Diagonal(x = diag(e2))))
986            }
987        }
988    }
989    
990    
991    ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"
992    ### ----------   the last 4: separately here
993    for(cl in diCls) {
994    setMethod("any", cl,
995              function (x, ..., na.rm) {
996                  if(any(x@Dim == 0)) FALSE
997                  else if(x@diag == "U") TRUE else any(x@x, ..., na.rm = na.rm)
998              })
999    setMethod("all",  cl, function (x, ..., na.rm) {
1000        n <- x@Dim[1]
1001        if(n >= 2) FALSE
1002        else if(n == 0 || x@diag == "U") TRUE
1003        else all(x@x, ..., na.rm = na.rm)
1004    })
1005    setMethod("prod", cl, function (x, ..., na.rm) {
1006        n <- x@Dim[1]
1007        if(n >= 2) 0
1008        else if(n == 0 || x@diag == "U") 1
1009        else ## n == 1, diag = "N" :
1010            prod(x@x, ..., na.rm = na.rm)
1011    })
1012    
1013    setMethod("sum", cl,
1014              function(x, ..., na.rm) {
1015                  r <- sum(x@x, ..., na.rm = na.rm)# double or integer, correctly
1016                  if(x@diag == "U" && !is.na(r)) r + x@Dim[1] else r
1017              })
1018    }
1019    
1020    ## The remaining ones are  max, min, range :
1021    
1022    setMethod("Summary", "ddiMatrix",
1023              function(x, ..., na.rm) {
1024                  if(any(x@Dim == 0)) callGeneric(numeric(0), ..., na.rm=na.rm)
1025                  else if(x@diag == "U")
1026                      callGeneric(x@x, 0, 1, ..., na.rm=na.rm)
1027                  else callGeneric(x@x, 0, ..., na.rm=na.rm)
1028              })
1029    setMethod("Summary", "ldiMatrix",
1030              function(x, ..., na.rm) {
1031                  if(any(x@Dim == 0)) callGeneric(logical(0), ..., na.rm=na.rm)
1032                  else if(x@diag == "U")
1033                      callGeneric(x@x, FALSE, TRUE, ..., na.rm=na.rm)
1034                  else callGeneric(x@x, FALSE, ..., na.rm=na.rm)
1035              })
1036    
1037    
1038    
# Line 669  Line 1047 
1047      invisible(x)      invisible(x)
1048  }  }
1049    
1050    ## somewhat consistent with "print" for sparseMatrix :
1051    setMethod("print", signature(x = "diagonalMatrix"), prDiag)
1052    
1053  setMethod("show", signature(object = "diagonalMatrix"),  setMethod("show", signature(object = "diagonalMatrix"),
1054            function(object) {            function(object) {
1055                d <- dim(object)                d <- dim(object)
1056                cl <- class(object)                cl <- class(object)
1057                cat(sprintf('%d x %d diagonal matrix of class "%s"\n',                cat(sprintf('%d x %d diagonal matrix of class "%s"',
1058                            d[1], d[2], cl))                            d[1], d[2], cl))
1059                  if(d[1] < 50) {
1060                      cat("\n")
1061                prDiag(object)                prDiag(object)
1062                  } else {
1063                      cat(", with diagonal entries\n")
1064                      show(diag(object))
1065                      invisible(object)
1066                  }
1067            })            })

Legend:
Removed from v.2157  
changed lines
  Added in v.2608

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