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 2904, Tue Sep 10 19:43:53 2013 UTC revision 2998, Tue Aug 12 10:27:15 2014 UTC
# Line 227  Line 227 
227                  n))                  n))
228  }  }
229    
230  ## diagonal -> triangular,  upper / lower depending on "partner":  ## diagonal -> triangular,  upper / lower depending on "partner" 'x':
231  diag2tT.u <- function(d, x, kind = .M.kind(d))  diag2tT.u <- function(d, x, kind = .M.kind(d))
232      .diag2tT(d, uplo = if(is(x,"triangularMatrix")) x@uplo else "U", kind)      .diag2tT(d, uplo = if(is(x,"triangularMatrix")) x@uplo else "U", kind)
233    
# Line 251  Line 251 
251  ## "hack"  instead of signature  x = "diagonalMatrix" :  ## "hack"  instead of signature  x = "diagonalMatrix" :
252  ##  ##
253  ## ddi*:  ## ddi*:
254  diag2tT <- function(from) .diag2tT(from, "U", "d")  di2tT <- function(from) .diag2tT(from, "U", "d")
255  setAs("ddiMatrix", "triangularMatrix", diag2tT)  setAs("ddiMatrix", "triangularMatrix", di2tT)
256  ##_no_longer_ setAs("ddiMatrix", "sparseMatrix", diag2tT)  ##_no_longer_ setAs("ddiMatrix", "sparseMatrix", di2tT)
257  ## needed too (otherwise <dense> -> Tsparse is taken):  ## needed too (otherwise <dense> -> Tsparse is taken):
258  setAs("ddiMatrix", "TsparseMatrix", diag2tT)  setAs("ddiMatrix", "TsparseMatrix", di2tT)
259  setAs("ddiMatrix", "dsparseMatrix", diag2tT)  setAs("ddiMatrix", "dsparseMatrix", di2tT)
260  setAs("ddiMatrix", "CsparseMatrix",  setAs("ddiMatrix", "CsparseMatrix",
261        function(from) as(.diag2tT(from, "U", "d"), "CsparseMatrix"))        function(from) as(.diag2tT(from, "U", "d"), "CsparseMatrix"))
262  setAs("ddiMatrix", "symmetricMatrix",  setAs("ddiMatrix", "symmetricMatrix", function(from) .diag2sT(from, "U", "d"))
       function(from) .diag2sT(from, "U", "d"))  
263  ##  ##
264  ## ldi*:  ## ldi*:
265  diag2tT <- function(from) .diag2tT(from, "U", "l")  di2tT <- function(from) .diag2tT(from, "U", "l")
266  setAs("ldiMatrix", "triangularMatrix", diag2tT)  setAs("ldiMatrix", "triangularMatrix", di2tT)
267  ##_no_longer_ setAs("ldiMatrix", "sparseMatrix", diag2tT)  ##_no_longer_ setAs("ldiMatrix", "sparseMatrix", di2tT)
268  ## needed too (otherwise <dense> -> Tsparse is taken):  ## needed too (otherwise <dense> -> Tsparse is taken):
269  setAs("ldiMatrix", "TsparseMatrix", diag2tT)  setAs("ldiMatrix", "TsparseMatrix", di2tT)
270  setAs("ldiMatrix", "lsparseMatrix", diag2tT)  setAs("ldiMatrix", "lsparseMatrix", di2tT)
271  setAs("ldiMatrix", "CsparseMatrix",  setAs("ldiMatrix", "CsparseMatrix",
272        function(from) as(.diag2tT(from, "U", "l"), "CsparseMatrix"))        function(from) as(.diag2tT(from, "U", "l"), "CsparseMatrix"))
273  setAs("ldiMatrix", "symmetricMatrix",  setAs("ldiMatrix", "symmetricMatrix", function(from) .diag2sT(from, "U", "l"))
274        function(from) .diag2sT(from, "U", "l"))  rm(di2tT)
   
275    
276  setAs("diagonalMatrix", "nMatrix",  setAs("diagonalMatrix", "nMatrix",
277        function(from) {        function(from) {
# Line 285  Line 283 
283    
284  setAs("diagonalMatrix", "nsparseMatrix", function(from) as(from, "nMatrix"))  setAs("diagonalMatrix", "nsparseMatrix", function(from) as(from, "nMatrix"))
285    
286  ## 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"
287  mkDiag <- function(x, n) {  mkDiag <- function(x, n) {
288      y <- matrix(as0(mod=mode(x)), n,n)      y <- matrix(as0(mod=mode(x)), n,n)
289      if (n > 0) y[1L + 0:(n - 1L) * (n + 1)] <- x      if (n > 0) y[1L + 0:(n - 1L) * (n + 1)] <- x
290      y      y
291  }  }
292    ## NB: diag(x,n) is really faster for n >= 20, and even more for large n
293    ## --> using diag() where possible, ==> .ddi2mat()
294    
295  setAs("diagonalMatrix", "matrix",  .diag2mat <- function(from)
296        function(from) {      ## want "ldiMatrix" -> <logical> "matrix"  (but integer -> <double> for now)
297            ## want "ldiMatrix" -> <logical> "matrix" :      mkDiag(if(from@diag == "U") as1(from@x) else from@x, n = from@Dim[1])
298            mkDiag(if(from@diag == "U") as1(from@x) else from@x,  
299                   n = from@Dim[1])  .ddi2mat <- function(from)
300        })      base::diag(if(from@diag == "U") as1(from@x) else from@x, nrow = from@Dim[1])
301    
302    setAs("ddiMatrix", "matrix", .ddi2mat)
303    ## the non-ddi diagonalMatrix -- only "ldiMatrix" currently:
304    setAs("diagonalMatrix", "matrix", .diag2mat)
305    
306  setMethod("as.vector", signature(x = "diagonalMatrix", mode="missing"),  setMethod("as.vector", signature(x = "diagonalMatrix", mode="missing"),
307            function(x, mode) {            function(x, mode) {
# Line 581  Line 585 
585            diagdiagprod, valueClass = "ddiMatrix")            diagdiagprod, valueClass = "ddiMatrix")
586    
587    
588    ## analogous to matdiagprod() below:
589  diagmatprod <- function(x, y) {  diagmatprod <- function(x, y) {
590      ## x is diagonalMatrix      ## x is diagonalMatrix
591      dx <- dim(x)      if(x@Dim[2] != nrow(y)) stop("non-matching dimensions")
592      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")  
593  }  }
594  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),
595            diagmatprod)            diagmatprod)
596  ## sneaky .. :  ## sneaky .. :
597  formals(diagmatprod) <- alist(x=, y=NULL)  formals(diagmatprod) <- alist(x=, y=NULL)
598  setMethod("crossprod", signature(x = "diagonalMatrix", y = "matrix"),  setMethod("crossprod",  signature(x = "diagonalMatrix", y = "matrix"), diagmatprod)
           diagmatprod)  
599    
600  diagGeprod <- function(x, y) {  diagGeprod <- function(x, y) {
601      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")  
602      if(x@diag != "U")      if(x@diag != "U")
603          y@x <- x@x * y@x          y@x <- x@x * y@x
604      y      y
# Line 611  Line 611 
611  setMethod("crossprod", signature(x = "diagonalMatrix", y = "lgeMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "lgeMatrix"),
612            diagGeprod)            diagGeprod)
613    
614    ## analogous to diagmatprod() above:
615  matdiagprod <- function(x, y) {  matdiagprod <- function(x, y) {
616      dx <- dim(x)      dx <- dim(x)
617      dy <- dim(y)      if(dx[2] != y@Dim[1]) stop("non-matching dimensions")
     if(dx[2] != dy[1]) stop("non-matching dimensions")  
618      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]))
619  }  }
620  setMethod("%*%", signature(x = "matrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "matrix", y = "diagonalMatrix"), matdiagprod)
           matdiagprod)  
621  formals(matdiagprod) <- alist(x=, y=NULL)  formals(matdiagprod) <- alist(x=, y=NULL)
622  setMethod("tcrossprod", signature(x = "matrix", y = "diagonalMatrix"),  setMethod("tcrossprod",signature(x = "matrix", y = "diagonalMatrix"), matdiagprod)
623            matdiagprod)  setMethod("crossprod", signature(x = "matrix", y = "diagonalMatrix"),
624              function(x, y=NULL) {
625                  dx <- dim(x)
626                  if(dx[1] != y@Dim[1]) stop("non-matching dimensions")
627                  Matrix(t(rep.int(y@x, dx[2]) * x))
628              })
629    
630  gediagprod <- function(x, y) {  gediagprod <- function(x, y) {
631      dx <- dim(x)      dx <- dim(x)
632      dy <- dim(y)      if(dx[2] != y@Dim[1]) stop("non-matching dimensions")
     if(dx[2] != dy[1]) stop("non-matching dimensions")  
633      if(y@diag == "N")      if(y@diag == "N")
634          x@x <- x@x * rep(y@x, each = dx[1])          x@x <- x@x * rep(y@x, each = dx[1])
635      x      x
# Line 654  Line 657 
657  ##        function(x, y = NULL) {  ##        function(x, y = NULL) {
658  ##           })  ##           })
659    
660    ##' @param x CsparseMatrix
661    ##' @param y diagonalMatrix
662    ##' @return x %*% y
663  Cspdiagprod <- function(x, y) {  Cspdiagprod <- function(x, y) {
664      dx <- dim(x <- .Call(Csparse_diagU2N, x))      m <- ncol(x <- .Call(Csparse_diagU2N, x))
665      dy <- dim(y)      if(m != y@Dim[1]) stop("non-matching dimensions")
666      if(dx[2] != dy[1]) stop("non-matching dimensions")      if(y@diag == "N") { ## otherwise: y == Diagonal(n) : multiplication is identity
     if(y@diag == "N") {  
667          if(!all(y@x[1L] == y@x[-1L]) && is(x, "symmetricMatrix"))          if(!all(y@x[1L] == y@x[-1L]) && is(x, "symmetricMatrix"))
668              x <- as(x, "generalMatrix")              x <- as(x, "generalMatrix")
669          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])
670          x@x <- x@x * y@x[ind]          x@x <- x@x * y@x[ind]
671      }          if(.hasSlot(x, "factors") && length(x@factors)) {# drop cashed ones
     if(is(x, "compMatrix") && length(xf <- x@factors)) {  
672          ## instead of dropping all factors, be smart about some          ## instead of dropping all factors, be smart about some
673          ## TODO ......          ## TODO ......
674          x@factors <- list()          x@factors <- list()
675      }      }
676        }
677      x      x
678  }  }
679    
680    ##' @param x diagonalMatrix
681    ##' @param y CsparseMatrix
682    ##' @return x %*% y
683  diagCspprod <- function(x, y) {  diagCspprod <- function(x, y) {
684      dx <- dim(x)      y <- .Call(Csparse_diagU2N, y)
685      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")  
686      if(x@diag == "N") {      if(x@diag == "N") {
687          if(!all(x@x[1L] == x@x[-1L]) && is(y, "symmetricMatrix"))          if(!all(x@x[1L] == x@x[-1L]) && is(y, "symmetricMatrix"))
688              y <- as(y, "generalMatrix")              y <- as(y, "generalMatrix")
689          y@x <- y@x * x@x[y@i + 1L]          y@x <- y@x * x@x[y@i + 1L]
690      }          if(.hasSlot(y, "factors") && length(yf <- y@factors)) {
     if(is(y, "compMatrix") && length(yf <- y@factors)) {  
         ## instead of dropping all factors, be smart about some  
691          ## TODO          ## TODO
692                if(FALSE) { ## instead of dropping all factors, be smart about some
693          keep <- character()          keep <- character()
694          if(iLU <- names(yf) == "LU") {                  if(any(iLU <- names(yf) == "LU")) {
695              ## TODO keep <- "LU"                      keep <- "LU"
696          }          }
697          y@factors <- yf[keep]          y@factors <- yf[keep]
698                } else y@factors <- list() ## for now
699            }
700      }      }
701      y      y
702  }  }
# Line 722  Line 730 
730  setMethod("%*%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
731            function(x, y) diagCspprod(x, y))            function(x, y) diagCspprod(x, y))
732    
733    ## instead of "sparseMatrix", use: [RT]sparse.. ("closer" in method dispatch)
734    for(cl in c("TsparseMatrix", "RsparseMatrix")) {
735    
736  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
737            function(x, y) diagCspprod(as(x, "CsparseMatrix"), y))            function(x, y) diagCspprod(as(x, "CsparseMatrix"), y))
738    
739  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
740            function(x, y) Cspdiagprod(as(x, "CsparseMatrix"), y))            function(x, y) Cspdiagprod(as(x, "CsparseMatrix"), y))
741    }
742    
743  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
744            function(x, y) Cspdiagprod(x, y))            function(x, y) Cspdiagprod(x, y))
# Line 907  Line 919 
919            function(e1,e2) {            function(e1,e2) {
920                n <- e1@Dim[1]                n <- e1@Dim[1]
921                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
922                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
923                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
924                    if(e1@diag == "U") {                    if(e1@diag == "U") {
925                        if(any((r <- callGeneric(1, e2)) != 1)) {                        if(any((r <- callGeneric(1, e2)) != 1)) {
# Line 929  Line 941 
941            function(e1,e2) {            function(e1,e2) {
942                n <- e2@Dim[1]                n <- e2@Dim[1]
943                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
944                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
945                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
946                    if(e2@diag == "U") {                    if(e2@diag == "U") {
947                        if(any((r <- callGeneric(e1, 1)) != 1)) {                        if(any((r <- callGeneric(e1, 1)) != 1)) {
# Line 952  Line 964 
964            function(e1,e2) {            function(e1,e2) {
965                n <- e1@Dim[1]                n <- e1@Dim[1]
966                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
967                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
968                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
969                    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)  
970                    ## storage.mode(E@x) <- "double"                    ## storage.mode(E@x) <- "double"
971                    if(e1@diag == "U") {                    if(e1@diag == "U") {
972                        if(any((r <- callGeneric(1, e2)) != 1)) {                        if(any((r <- callGeneric(1, e2)) != 1)) {
# Line 977  Line 988 
988            function(e1,e2) {            function(e1,e2) {
989                n <- e2@Dim[1]                n <- e2@Dim[1]
990                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
991                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
992                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
993                    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)  
994                    ## storage.mode(E@x) <- "double"                    ## storage.mode(E@x) <- "double"
995                    if(e2@diag == "U") {                    if(e2@diag == "U") {
996                        if(any((r <- callGeneric(e1, 1)) != 1)) {                        if(any((r <- callGeneric(e1, 1)) != 1)) {
# Line 1010  Line 1020 
1020            function(e1,e2) {            function(e1,e2) {
1021                n <- e1@Dim[1]                n <- e1@Dim[1]
1022                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
1023                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
1024                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
1025                    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)  
1026                    ## storage.mode(E@x) <- "logical"                    ## storage.mode(E@x) <- "logical"
1027                    if(e1@diag == "U") {                    if(e1@diag == "U") {
1028                        if(any((r <- callGeneric(1, e2)) != 1)) {                        if(any((r <- callGeneric(1, e2)) != 1)) {
# Line 1036  Line 1045 
1045            function(e1,e2) {            function(e1,e2) {
1046                n <- e1@Dim[1]                n <- e1@Dim[1]
1047                f0 <- callGeneric(FALSE, e2)                f0 <- callGeneric(FALSE, e2)
1048                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
1049                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
1050    
1051                    if(e1@diag == "U") {                    if(e1@diag == "U") {
# Line 1091  Line 1100 
1100      }      }
1101  }  }
1102    
   
1103  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"
1104  ### ----------   the last 4: separately here  ### ----------   the last 4: separately here
1105  for(cl in diCls) {  for(cl in diCls) {
# Line 1170  Line 1178 
1178                }                }
1179            })            })
1180    
1181  rm(arg1, arg2, other, DI, cl, c1, c2,  rm(arg1, arg2, other, DI, Fun, cl, c1, c2,
1182     dense.subCl, diCls)# not used elsewhere     dense.subCl, diCls)# not used elsewhere
1183    
1184  setMethod("summary", signature(object = "diagonalMatrix"),  setMethod("summary", signature(object = "diagonalMatrix"),

Legend:
Removed from v.2904  
changed lines
  Added in v.2998

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