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 3002, Mon Sep 1 15:31:40 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 658  Line 661 
661  ##' @param y diagonalMatrix  ##' @param y diagonalMatrix
662  ##' @return x %*% y  ##' @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")
     if(dx[2] != dy[1]) stop("non-matching dimensions")  
666      if(y@diag == "N") { ## otherwise: y == Diagonal(n) : multiplication is identity      if(y@diag == "N") { ## otherwise: y == Diagonal(n) : multiplication is identity
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(is(x, "compMatrix") && length(xf <- x@factors)) {          if(.hasSlot(x, "factors") && length(x@factors)) {# drop cashed ones
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()
# Line 679  Line 681 
681  ##' @param y CsparseMatrix  ##' @param y CsparseMatrix
682  ##' @return x %*% y  ##' @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(is(y, "compMatrix") && length(yf <- y@factors)) {          if(.hasSlot(y, "factors") && length(yf <- y@factors)) {
691              ## TODO              ## TODO
692              if(FALSE) { ## instead of dropping all factors, be smart about some              if(FALSE) { ## instead of dropping all factors, be smart about some
693                  keep <- character()                  keep <- character()
# Line 918  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 940  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 963  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 988  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 1021  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 1047  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 1102  Line 1100 
1100      }      }
1101  }  }
1102    
1103    ## Group methods "Math", "Math2" in                     --> ./Math.R
1104    
1105  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"
1106  ### ----------   the last 4: separately here  ### ----------   the last 4: separately here
# Line 1181  Line 1180 
1180                }                }
1181            })            })
1182    
1183  rm(arg1, arg2, other, DI, cl, c1, c2,  rm(arg1, arg2, other, DI, Fun, cl, c1, c2,
1184     dense.subCl, diCls)# not used elsewhere     dense.subCl, diCls)# not used elsewhere
1185    
1186  setMethod("summary", signature(object = "diagonalMatrix"),  setMethod("summary", signature(object = "diagonalMatrix"),

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

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