SCM

SCM Repository

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

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

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

revision 2912, Sat Sep 14 17:09:49 2013 UTC revision 3046, Mon Feb 16 14:40:37 2015 UTC
# Line 198  Line 198 
198      as(.bdiag(alis), "CsparseMatrix")      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)) {  .diag2tT <- function(from, uplo = "U", kind = .M.kind(from)) {
210      ## to triangular Tsparse      ## to triangular Tsparse
# Line 227  Line 234 
234                  n))                  n))
235  }  }
236    
237  ## diagonal -> triangular,  upper / lower depending on "partner":  ## diagonal -> triangular,  upper / lower depending on "partner" 'x':
238  diag2tT.u <- function(d, x, kind = .M.kind(d))  diag2tT.u <- function(d, x, kind = .M.kind(d))
239      .diag2tT(d, uplo = if(is(x,"triangularMatrix")) x@uplo else "U", kind)      .diag2tT(d, uplo = if(is(x,"triangularMatrix")) x@uplo else "U", kind)
240    
# Line 251  Line 258 
258  ## "hack"  instead of signature  x = "diagonalMatrix" :  ## "hack"  instead of signature  x = "diagonalMatrix" :
259  ##  ##
260  ## ddi*:  ## ddi*:
261  diag2tT <- function(from) .diag2tT(from, "U", "d")  di2tT <- function(from) .diag2tT(from, "U", "d")
262  setAs("ddiMatrix", "triangularMatrix", diag2tT)  setAs("ddiMatrix", "triangularMatrix", di2tT)
263  ##_no_longer_ setAs("ddiMatrix", "sparseMatrix", diag2tT)  ##_no_longer_ setAs("ddiMatrix", "sparseMatrix", di2tT)
264  ## needed too (otherwise <dense> -> Tsparse is taken):  ## needed too (otherwise <dense> -> Tsparse is taken):
265  setAs("ddiMatrix", "TsparseMatrix", diag2tT)  setAs("ddiMatrix", "TsparseMatrix", di2tT)
266  setAs("ddiMatrix", "dsparseMatrix", diag2tT)  setAs("ddiMatrix", "dsparseMatrix", di2tT)
267  setAs("ddiMatrix", "CsparseMatrix",  setAs("ddiMatrix", "CsparseMatrix",
268        function(from) as(.diag2tT(from, "U", "d"), "CsparseMatrix"))        function(from) as(.diag2tT(from, "U", "d"), "CsparseMatrix"))
269  setAs("ddiMatrix", "symmetricMatrix",  setAs("ddiMatrix", "symmetricMatrix", function(from) .diag2sT(from, "U", "d"))
       function(from) .diag2sT(from, "U", "d"))  
270  ##  ##
271  ## ldi*:  ## ldi*:
272  diag2tT <- function(from) .diag2tT(from, "U", "l")  di2tT <- function(from) .diag2tT(from, "U", "l")
273  setAs("ldiMatrix", "triangularMatrix", diag2tT)  setAs("ldiMatrix", "triangularMatrix", di2tT)
274  ##_no_longer_ setAs("ldiMatrix", "sparseMatrix", diag2tT)  ##_no_longer_ setAs("ldiMatrix", "sparseMatrix", di2tT)
275  ## needed too (otherwise <dense> -> Tsparse is taken):  ## needed too (otherwise <dense> -> Tsparse is taken):
276  setAs("ldiMatrix", "TsparseMatrix", diag2tT)  setAs("ldiMatrix", "TsparseMatrix", di2tT)
277  setAs("ldiMatrix", "lsparseMatrix", diag2tT)  setAs("ldiMatrix", "lsparseMatrix", di2tT)
278  setAs("ldiMatrix", "CsparseMatrix",  setAs("ldiMatrix", "CsparseMatrix",
279        function(from) as(.diag2tT(from, "U", "l"), "CsparseMatrix"))        function(from) as(.diag2tT(from, "U", "l"), "CsparseMatrix"))
280  setAs("ldiMatrix", "symmetricMatrix",  setAs("ldiMatrix", "symmetricMatrix", function(from) .diag2sT(from, "U", "l"))
281        function(from) .diag2sT(from, "U", "l"))  rm(di2tT)
   
282    
283  setAs("diagonalMatrix", "nMatrix",  setAs("diagonalMatrix", "nMatrix",
284        function(from) {        function(from) {
           n <- from@Dim[1]  
285            i <- if(from@diag == "U") integer(0) else which(isN0(from@x)) - 1L            i <- if(from@diag == "U") integer(0) else which(isN0(from@x)) - 1L
286            new("ntTMatrix", i = i, j = i, diag = from@diag,            new("ntTMatrix", i = i, j = i, diag = from@diag,
287                Dim = from@Dim, Dimnames = from@Dimnames)                Dim = from@Dim, Dimnames = from@Dimnames)
# Line 285  Line 289 
289    
290  setAs("diagonalMatrix", "nsparseMatrix", function(from) as(from, "nMatrix"))  setAs("diagonalMatrix", "nsparseMatrix", function(from) as(from, "nMatrix"))
291    
292  ## Cheap fast substitute for diag() which *does* preserve the mode of x :  ##' A version of diag(x,n) which *does* preserve the mode of x, where diag() "fails"
293  mkDiag <- function(x, n) {  mkDiag <- function(x, n) {
294      y <- matrix(as0(mod=mode(x)), n,n)      y <- matrix(as0(mod=mode(x)), n,n)
295      if (n > 0) y[1L + 0:(n - 1L) * (n + 1)] <- x      if (n > 0) y[1L + 0:(n - 1L) * (n + 1)] <- x
296      y      y
297  }  }
298    ## NB: diag(x,n) is really faster for n >= 20, and even more for large n
299    ## --> using diag() where possible, ==> .ddi2mat()
300    
301  setAs("diagonalMatrix", "matrix",  .diag2mat <- function(from)
302        function(from) {      ## want "ldiMatrix" -> <logical> "matrix"  (but integer -> <double> for now)
303            ## want "ldiMatrix" -> <logical> "matrix" :      mkDiag(if(from@diag == "U") as1(from@x) else from@x, n = from@Dim[1])
304            mkDiag(if(from@diag == "U") as1(from@x) else from@x,  
305                   n = from@Dim[1])  .ddi2mat <- function(from)
306        })      base::diag(if(from@diag == "U") as1(from@x) else from@x, nrow = from@Dim[1])
307    
308    setAs("ddiMatrix", "matrix", .ddi2mat)
309    ## the non-ddi diagonalMatrix -- only "ldiMatrix" currently:
310    setAs("diagonalMatrix", "matrix", .diag2mat)
311    
312  setMethod("as.vector", signature(x = "diagonalMatrix", mode="missing"),  setMethod("as.vector", signature(x = "diagonalMatrix", mode="missing"),
313            function(x, mode) {            function(x, mode) {
# Line 516  Line 526 
526            function(x) { x@Dimnames <- x@Dimnames[2:1] ; x })            function(x) { x@Dimnames <- x@Dimnames[2:1] ; x })
527    
528  setMethod("isDiagonal",   "diagonalMatrix", function(object) TRUE)  setMethod("isDiagonal",   "diagonalMatrix", function(object) TRUE)
529  setMethod("isTriangular", "diagonalMatrix", function(object, ...) TRUE)  setMethod("isTriangular", "diagonalMatrix", function(object, upper=NA, ...) TRUE)
530  setMethod("isSymmetric",  "diagonalMatrix", function(object, ...) TRUE)  setMethod("isSymmetric",  "diagonalMatrix", function(object, ...) TRUE)
531    
532  setMethod("symmpart", signature(x = "diagonalMatrix"), function(x) x)  setMethod("symmpart", signature(x = "diagonalMatrix"), function(x) x)
533  setMethod("skewpart", signature(x = "diagonalMatrix"), setZero)  setMethod("skewpart", signature(x = "diagonalMatrix"), function(x) .setZero(x))
534    
535  setMethod("chol", signature(x = "ddiMatrix"),  setMethod("chol", signature(x = "ddiMatrix"),
536            function(x, pivot, ...) {            function(x, pivot, ...) {
# Line 581  Line 591 
591            diagdiagprod, valueClass = "ddiMatrix")            diagdiagprod, valueClass = "ddiMatrix")
592    
593    
594    ## analogous to matdiagprod() below:
595  diagmatprod <- function(x, y) {  diagmatprod <- function(x, y) {
596      ## x is diagonalMatrix      ## x is diagonalMatrix
597      dx <- dim(x)      if(x@Dim[2] != nrow(y)) stop("non-matching dimensions")
598      dy <- dim(y)      Matrix(if(x@diag == "U") y else x@x * y)
     if(dx[2] != dy[1]) stop("non-matching dimensions")  
     as(if(x@diag == "U") y else x@x * y, "Matrix")  
599  }  }
600  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),
601            diagmatprod)            diagmatprod)
602  ## sneaky .. :  ## sneaky .. :
603  formals(diagmatprod) <- alist(x=, y=NULL)  formals(diagmatprod) <- alist(x=, y=NULL)
604  setMethod("crossprod", signature(x = "diagonalMatrix", y = "matrix"),  setMethod("crossprod",  signature(x = "diagonalMatrix", y = "matrix"), diagmatprod)
           diagmatprod)  
605    
606  diagGeprod <- function(x, y) {  diagGeprod <- function(x, y) {
607      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")  
608      if(x@diag != "U")      if(x@diag != "U")
609          y@x <- x@x * y@x          y@x <- x@x * y@x
610      y      y
# Line 611  Line 617 
617  setMethod("crossprod", signature(x = "diagonalMatrix", y = "lgeMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "lgeMatrix"),
618            diagGeprod)            diagGeprod)
619    
620    ## analogous to diagmatprod() above:
621  matdiagprod <- function(x, y) {  matdiagprod <- function(x, y) {
622      dx <- dim(x)      dx <- dim(x)
623      dy <- dim(y)      if(dx[2] != y@Dim[1]) stop("non-matching dimensions")
     if(dx[2] != dy[1]) stop("non-matching dimensions")  
624      Matrix(if(y@diag == "U") x else x * rep(y@x, each = dx[1]))      Matrix(if(y@diag == "U") x else x * rep(y@x, each = dx[1]))
625  }  }
626  setMethod("%*%", signature(x = "matrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "matrix", y = "diagonalMatrix"), matdiagprod)
           matdiagprod)  
627  formals(matdiagprod) <- alist(x=, y=NULL)  formals(matdiagprod) <- alist(x=, y=NULL)
628  setMethod("tcrossprod", signature(x = "matrix", y = "diagonalMatrix"),  setMethod("tcrossprod",signature(x = "matrix", y = "diagonalMatrix"), matdiagprod)
629            matdiagprod)  setMethod("crossprod", signature(x = "matrix", y = "diagonalMatrix"),
630              function(x, y=NULL) {
631                  dx <- dim(x)
632                  if(dx[1] != y@Dim[1]) stop("non-matching dimensions")
633                  Matrix(t(rep.int(y@x, dx[2]) * x))
634              })
635    
636  gediagprod <- function(x, y) {  gediagprod <- function(x, y) {
637      dx <- dim(x)      dx <- dim(x)
638      dy <- dim(y)      if(dx[2] != y@Dim[1]) stop("non-matching dimensions")
     if(dx[2] != dy[1]) stop("non-matching dimensions")  
639      if(y@diag == "N")      if(y@diag == "N")
640          x@x <- x@x * rep(y@x, each = dx[1])          x@x <- x@x * rep(y@x, each = dx[1])
641      x      x
# Line 658  Line 667 
667  ##' @param y diagonalMatrix  ##' @param y diagonalMatrix
668  ##' @return x %*% y  ##' @return x %*% y
669  Cspdiagprod <- function(x, y) {  Cspdiagprod <- function(x, y) {
670      dx <- dim(x <- .Call(Csparse_diagU2N, x))      m <- ncol(x <- .Call(Csparse_diagU2N, x))
671      dy <- dim(y)      if(m != y@Dim[1]) stop("non-matching dimensions")
     if(dx[2] != dy[1]) stop("non-matching dimensions")  
672      if(y@diag == "N") { ## otherwise: y == Diagonal(n) : multiplication is identity      if(y@diag == "N") { ## otherwise: y == Diagonal(n) : multiplication is identity
673          if(!all(y@x[1L] == y@x[-1L]) && is(x, "symmetricMatrix"))          if(!all(y@x[1L] == y@x[-1L]) && is(x, "symmetricMatrix"))
674              x <- as(x, "generalMatrix")              x <- as(x, "generalMatrix")
675          ind <- rep.int(seq_len(dx[2]), x@p[-1] - x@p[-dx[2]-1L])          ind <- rep.int(seq_len(m), x@p[-1] - x@p[-m-1L])
676          x@x <- x@x * y@x[ind]          x@x <- x@x * y@x[ind]
677          if(is(x, "compMatrix") && length(xf <- x@factors)) {          if(.hasSlot(x, "factors") && length(x@factors)) {# drop cashed ones
678              ## instead of dropping all factors, be smart about some              ## instead of dropping all factors, be smart about some
679              ## TODO ......              ## TODO ......
680              x@factors <- list()              x@factors <- list()
# Line 679  Line 687 
687  ##' @param y CsparseMatrix  ##' @param y CsparseMatrix
688  ##' @return x %*% y  ##' @return x %*% y
689  diagCspprod <- function(x, y) {  diagCspprod <- function(x, y) {
690      dx <- dim(x)      y <- .Call(Csparse_diagU2N, y)
691      dy <- dim(y <- .Call(Csparse_diagU2N, y))      if(x@Dim[2] != y@Dim[1]) stop("non-matching dimensions")
     if(dx[2] != dy[1]) stop("non-matching dimensions")  
692      if(x@diag == "N") {      if(x@diag == "N") {
693          if(!all(x@x[1L] == x@x[-1L]) && is(y, "symmetricMatrix"))          if(!all(x@x[1L] == x@x[-1L]) && is(y, "symmetricMatrix"))
694              y <- as(y, "generalMatrix")              y <- as(y, "generalMatrix")
695          y@x <- y@x * x@x[y@i + 1L]          y@x <- y@x * x@x[y@i + 1L]
696          if(is(y, "compMatrix") && length(yf <- y@factors)) {          if(.hasSlot(y, "factors") && length(y@factors)) {
697              ## TODO       ## if(.hasSlot(y, "factors") && length(yf <- y@factors)) { ## -- TODO? --
698              if(FALSE) { ## instead of dropping all factors, be smart about some              ## instead of dropping all factors, be smart about some
699                  keep <- character()              ## keep <- character()
700                  if(any(iLU <- names(yf) == "LU")) {              ## if(any(names(yf) == "LU")) { ## <- not easy: y = P'LUQ,  x y = xP'LUQ => LU ???
701                      keep <- "LU"              ##     keep <- "LU"
702                  }              ## }
703                  y@factors <- yf[keep]              ## y@factors <- yf[keep]
704              } else y@factors <- list() ## for now              y@factors <- list()
705          }          }
706      }      }
707      y      y
# Line 833  Line 840 
840      ## result must be triangular      ## result must be triangular
841      r <- callGeneric(d1 <- .diag.x(e1), diag(e2)) # error if not "compatible"      r <- callGeneric(d1 <- .diag.x(e1), diag(e2)) # error if not "compatible"
842      ## Check what happens with non-diagonals, i.e. (0 o 0), (FALSE o 0), ...:      ## Check what happens with non-diagonals, i.e. (0 o 0), (FALSE o 0), ...:
843      e1.0 <- if(.n1 <- is.numeric(d1   )) 0 else FALSE      e1.0 <- if(is.numeric(d1)) 0 else FALSE
844      r00 <- callGeneric(e1.0, if(.n2 <- is.numeric(e2[0])) 0 else FALSE)      r00 <- callGeneric(e1.0, if(.n2 <- is.numeric(e2[0])) 0 else FALSE)
845      if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* triangular      if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* triangular
846          diag(e2) <- r          diag(e2) <- r
# Line 918  Line 925 
925            function(e1,e2) {            function(e1,e2) {
926                n <- e1@Dim[1]                n <- e1@Dim[1]
927                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
928                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
929                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
930                    if(e1@diag == "U") {                    if(e1@diag == "U") {
931                        if(any((r <- callGeneric(1, e2)) != 1)) {                        if(any((r <- callGeneric(1, e2)) != 1)) {
# Line 940  Line 947 
947            function(e1,e2) {            function(e1,e2) {
948                n <- e2@Dim[1]                n <- e2@Dim[1]
949                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
950                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
951                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
952                    if(e2@diag == "U") {                    if(e2@diag == "U") {
953                        if(any((r <- callGeneric(e1, 1)) != 1)) {                        if(any((r <- callGeneric(e1, 1)) != 1)) {
# Line 963  Line 970 
970            function(e1,e2) {            function(e1,e2) {
971                n <- e1@Dim[1]                n <- e1@Dim[1]
972                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
973                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
974                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
975                    E <- copyClass(e1, "ddiMatrix", c("diag", "Dim", "Dimnames"))#FIXME: if ok, check=FALSE                    E <- copyClass(e1, "ddiMatrix", c("diag", "Dim", "Dimnames"), check=FALSE)
                   ## E <- copyClass(e1, "ddiMatrix", check=FALSE)  
976                    ## storage.mode(E@x) <- "double"                    ## storage.mode(E@x) <- "double"
977                    if(e1@diag == "U") {                    if(e1@diag == "U") {
978                        if(any((r <- callGeneric(1, e2)) != 1)) {                        if(any((r <- callGeneric(1, e2)) != 1)) {
# Line 988  Line 994 
994            function(e1,e2) {            function(e1,e2) {
995                n <- e2@Dim[1]                n <- e2@Dim[1]
996                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
997                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
998                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
999                    E <- copyClass(e2, "ddiMatrix", c("diag", "Dim", "Dimnames"))#FIXME: if ok, check=FALSE                    E <- copyClass(e2, "ddiMatrix", c("diag", "Dim", "Dimnames"), check=FALSE)
                   ## E <- copyClass(e2, "ddiMatrix", check=FALSE)  
1000                    ## storage.mode(E@x) <- "double"                    ## storage.mode(E@x) <- "double"
1001                    if(e2@diag == "U") {                    if(e2@diag == "U") {
1002                        if(any((r <- callGeneric(e1, 1)) != 1)) {                        if(any((r <- callGeneric(e1, 1)) != 1)) {
# Line 1021  Line 1026 
1026            function(e1,e2) {            function(e1,e2) {
1027                n <- e1@Dim[1]                n <- e1@Dim[1]
1028                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
1029                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
1030                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
1031                    E <- copyClass(e1, "ldiMatrix", c("diag", "Dim", "Dimnames"))#FIXME: if ok, check=FALSE                    E <- copyClass(e1, "ldiMatrix", c("diag", "Dim", "Dimnames"), check=FALSE)
                   ## E <- copyClass(e1, "ldiMatrix", check=FALSE)  
1032                    ## storage.mode(E@x) <- "logical"                    ## storage.mode(E@x) <- "logical"
1033                    if(e1@diag == "U") {                    if(e1@diag == "U") {
1034                        if(any((r <- callGeneric(1, e2)) != 1)) {                        if(any((r <- callGeneric(1, e2)) != 1)) {
# Line 1047  Line 1051 
1051            function(e1,e2) {            function(e1,e2) {
1052                n <- e1@Dim[1]                n <- e1@Dim[1]
1053                f0 <- callGeneric(FALSE, e2)                f0 <- callGeneric(FALSE, e2)
1054                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
1055                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
1056    
1057                    if(e1@diag == "U") {                    if(e1@diag == "U") {
# Line 1102  Line 1106 
1106      }      }
1107  }  }
1108    
1109    ## Group methods "Math", "Math2" in                     --> ./Math.R
1110    
1111  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"
1112  ### ----------   the last 4: separately here  ### ----------   the last 4: separately here
# Line 1181  Line 1186 
1186                }                }
1187            })            })
1188    
1189  rm(arg1, arg2, other, DI, cl, c1, c2,  rm(arg1, arg2, other, DI, Fun, cl, c1, c2,
1190     dense.subCl, diCls)# not used elsewhere     dense.subCl, diCls)# not used elsewhere
1191    
1192  setMethod("summary", signature(object = "diagonalMatrix"),  setMethod("summary", signature(object = "diagonalMatrix"),

Legend:
Removed from v.2912  
changed lines
  Added in v.3046

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