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 2268, Mon Sep 22 15:17:12 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  .sparseDiagonal <- function(n, x = rep.int(1,n), uplo = "U", shape = "t") {  .sparseDiagonal <- function(n, x = rep.int(1,m), uplo = "U",
36                                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      stopifnot(is.character(shape), nchar(shape) == 1,      m <- length(cols)
43                any(shape == c("t","s","g"))) # triangular / symmetric / general      if(missing(kind))
44      kind <-      kind <-
45          if(is.double(x)) "d"          if(is.double(x)) "d"
46          else if(is.logical(x)) "l"          else if(is.logical(x)) "l"
# Line 45  Line 48 
48              storage.mode(x) <- "double"              storage.mode(x) <- "double"
49              "d"              "d"
50          }          }
51      new(paste(kind, shape, "CMatrix", sep=''),      else stopifnot(any(kind == c("d","l","n")))
52          Dim = c(n,n), x = x, uplo = uplo,      if(kind != "n") {
53          i = if(n) 0:(n - 1L) else integer(0), p = 0:n)          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()  ## Pkg 'spdep' had (relatively slow) versions of this as_dsCMatrix_I()
73  .symDiagonal <- function(n, x = rep.int(1,n), uplo = "U")  .symDiagonal <- function(n, x = rep.int(1,n), uplo = "U")
74      .sparseDiagonal(n, x, uplo, shape = "s")      .sparseDiagonal(n, x, uplo, shape = "s")
75    
76  ## instead of   diagU2N(as(Diagonal(n), "CsparseMatrix")), diag = "N" in any case:  # instead of   diagU2N(as(Diagonal(n), "CsparseMatrix")), diag = "N" in any case:
77  .trDiagonal <- function(n, x = rep.int(1,n), uplo = "U")  .trDiagonal <- function(n, x = rep.int(1,n), uplo = "U")
78      .sparseDiagonal(n, x, uplo, shape = "t")      .sparseDiagonal(n, x, uplo, shape = "t")
79    
# Line 217  Line 236 
236  ##_no_longer_ 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 228  Line 248 
248  ##_no_longer_ 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 295  Line 316 
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 319  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        })        })
# Line 334  Line 355 
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, "TsparseMatrix")      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 348  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 589  Line 616 
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, "TsparseMatrix"), 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, "TsparseMatrix")))            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, "TsparseMatrix"), 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, "TsparseMatrix")))            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, "TsparseMatrix") %*% y)            function(x, y) diagCspprod(as(x, "CsparseMatrix"), y))
683    
684  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
685            function(x, y) x %*% as(y, "TsparseMatrix"))            function(x, y) Cspdiagprod(as(x, "CsparseMatrix"), y))
686  ## 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)  
687  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
688            function(x, y) x %*% as(y, "CsparseMatrix"))            function(x, y) Cspdiagprod(x, y))
689  ## NB: this is *not* needed for Tsparse & Rsparse  
690  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*
691  ##       do indeed work by going through sparse (and *not* ddense)!  ##       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, ...) {
695                a@x <- 1/ a@x                a@x <- 1/ a@x
# Line 673  Line 739 
739          .diag.2N(e1)          .diag.2N(e1)
740      }      }
741      else { ## result not diagonal, but at least symmetric:      else { ## result not diagonal, but at least symmetric:
742            ## e.g., m == m
743          isNum <- (is.numeric(r) || is.numeric(r00))          isNum <- (is.numeric(r) || is.numeric(r00))
744          isLog <- (is.logical(r) || is.logical(r00))          isLog <- (is.logical(r) || is.logical(r00))
745            Matrix.msg("exploding <diag> o <diag> into dense matrix", .M.level = 2)
         if(getOption("verbose"))  
             message("exploding  <diag>  o  <diag>  into dense matrix")  
746          d <- e1@Dim          d <- e1@Dim
747          n <- d[1]          n <- d[1]
748          stopifnot(length(r) == n)          stopifnot(length(r) == n)
# Line 705  Line 770 
770              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)
771  }  }
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 891  Line 956 
956  for(other in c("ANY", "Matrix", "dMatrix")) {  for(other in c("ANY", "Matrix", "dMatrix")) {
957      ## ddi*:      ## ddi*:
958      setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other),      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 = other, 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 = other),      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 = other, 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) {  for(DI in diCls) {
975      for(c2 in c("denseMatrix", "Matrix")) {      for(c2 in c(dense.subCl, "Matrix")) {
976          for(Fun in c("*", "^", "&")) {          for(Fun in c("*", "^", "&")) {
977              setMethod(Fun, signature(e1 = DI, e2 = c2),              setMethod(Fun, signature(e1 = DI, e2 = c2),
978                        function(e1,e2) callGeneric(e1, Diagonal(x = diag(e2))))                        function(e1,e2) callGeneric(e1, Diagonal(x = diag(e2))))
# Line 918  Line 989 
989    
990    
991  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"
992  ### ----------  any, all: separately here  ### ----------   the last 4: separately here
993  for(cl in diCls) {  for(cl in diCls) {
994  setMethod("any", cl,  setMethod("any", cl,
995            function (x, ..., na.rm) {            function (x, ..., na.rm) {
996                if(any(x@Dim == 0)) FALSE                if(any(x@Dim == 0)) FALSE
997                else if(x@diag == "U") TRUE else any(x@x, ..., na.rm = na.rm)                else if(x@diag == "U") TRUE else any(x@x, ..., na.rm = na.rm)
998            })            })
999  setMethod("all",  cl, function (x, ..., na.rm) any(x@Dim == 0))  setMethod("all",  cl, function (x, ..., na.rm) {
1000  setMethod("prod", cl, function (x, ..., na.rm) as.numeric(any(x@Dim == 0)))      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,  setMethod("sum", cl,
1014            function(x, ..., na.rm) {            function(x, ..., na.rm) {

Legend:
Removed from v.2268  
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