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 2237, Fri Jul 25 06:55:42 2008 UTC revision 2363, Thu Apr 9 20:45:32 2009 UTC
# Line 214  Line 214 
214  ## ddi*:  ## ddi*:
215  diag2tT <- function(from) .diag2tT(from, "U", "d")  diag2tT <- function(from) .diag2tT(from, "U", "d")
216  setAs("ddiMatrix", "triangularMatrix", diag2tT)  setAs("ddiMatrix", "triangularMatrix", diag2tT)
217  setAs("ddiMatrix", "sparseMatrix", diag2tT)  ##_no_longer_ setAs("ddiMatrix", "sparseMatrix", diag2tT)
218  ## needed too (otherwise <dense> -> Tsparse is taken):  ## needed too (otherwise <dense> -> Tsparse is taken):
219  setAs("ddiMatrix", "TsparseMatrix", diag2tT)  setAs("ddiMatrix", "TsparseMatrix", diag2tT)
220  setAs("ddiMatrix", "CsparseMatrix",  setAs("ddiMatrix", "CsparseMatrix",
# Line 225  Line 225 
225  ## ldi*:  ## ldi*:
226  diag2tT <- function(from) .diag2tT(from, "U", "l")  diag2tT <- function(from) .diag2tT(from, "U", "l")
227  setAs("ldiMatrix", "triangularMatrix", diag2tT)  setAs("ldiMatrix", "triangularMatrix", diag2tT)
228  setAs("ldiMatrix", "sparseMatrix", diag2tT)  ##_no_longer_ setAs("ldiMatrix", "sparseMatrix", diag2tT)
229  ## needed too (otherwise <dense> -> Tsparse is taken):  ## needed too (otherwise <dense> -> Tsparse is taken):
230  setAs("ldiMatrix", "TsparseMatrix", diag2tT)  setAs("ldiMatrix", "TsparseMatrix", diag2tT)
231  setAs("ldiMatrix", "CsparseMatrix",  setAs("ldiMatrix", "CsparseMatrix",
# Line 272  Line 272 
272  setAs("diagonalMatrix", "generalMatrix", # prefer sparse:  setAs("diagonalMatrix", "generalMatrix", # prefer sparse:
273        function(from) as(as(from, "CsparseMatrix"), "generalMatrix"))        function(from) as(as(from, "CsparseMatrix"), "generalMatrix"))
274    
275    setAs("diagonalMatrix", "denseMatrix",
276          function(from) as(as(from, "CsparseMatrix"), "denseMatrix"))
277    
278  .diag.x <- function(m) if(m@diag == "U") rep.int(as1(m@x), m@Dim[1]) else m@x  .diag.x <- function(m) if(m@diag == "U") rep.int(as1(m@x), m@Dim[1]) else m@x
279    
280  .diag.2N <- function(m) {  .diag.2N <- function(m) {
# Line 292  Line 295 
295            d <- dim(from)            d <- dim(from)
296            if(d[1] != (n <- d[2])) stop("non-square matrix")            if(d[1] != (n <- d[2])) stop("non-square matrix")
297            if(any(from[row(from) != col(from)] != 0))            if(any(from[row(from) != col(from)] != 0))
298                stop("has non-zero off-diagonal entries")                stop("matrix with non-zero off-diagonals cannot be coerced to diagonalMatrix")
299            x <- diag(from)            x <- diag(from)
300            if(is.logical(x)) {            if(is.logical(x)) {
301                cl <- "ldiMatrix"                cl <- "ldiMatrix"
302                uni <- all(x)                uni <- allTrue(x) ## uni := {is it unit-diagonal ?}
303            } else {            } else {
304                cl <- "ddiMatrix"                cl <- "ddiMatrix"
305                uni <- all(x == 1)                uni <- allTrue(x == 1)
306                storage.mode(x) <- "double"                storage.mode(x) <- "double"
307            } ## TODO: complex            } ## TODO: complex
308            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 316  Line 319 
319            x <- diag(from)            x <- diag(from)
320            if(is.logical(x)) {            if(is.logical(x)) {
321                cl <- "ldiMatrix"                cl <- "ldiMatrix"
322                uni <- all(x)                uni <- allTrue(x)
323            } else {            } else {
324                cl <- "ddiMatrix"                cl <- "ddiMatrix"
325                uni <- all(x == 1)                uni <- allTrue(x == 1)
326                storage.mode(x) <- "double"                storage.mode(x) <- "double"
327            }            } ## TODO: complex
328            new(cl, Dim = c(n,n), diag = if(uni) "U" else "N",            new(cl, Dim = c(n,n), diag = if(uni) "U" else "N",
329                x = if(uni) x[FALSE] else x)                x = if(uni) x[FALSE] else x)
330        })        })
# Line 331  Line 334 
334            function(x = 1, nrow, ncol) .diag.x(x))            function(x = 1, nrow, ncol) .diag.x(x))
335    
336  subDiag <- function(x, i, j, ..., drop) {  subDiag <- function(x, i, j, ..., drop) {
337      x <- as(x, "sparseMatrix")      x <- as(x, "TsparseMatrix")
338      x <- if(missing(i))      x <- if(missing(i))
339          x[, j, drop=drop]          x[, j, drop=drop]
340      else if(missing(j))      else if(missing(j))
341          x[i, , drop=drop]          if(nargs() == 4) x[i, , drop=drop] else x[i, drop=drop]
342      else      else
343          x[i,j, drop=drop]          x[i,j, drop=drop]
344      if(isS4(x) && isDiagonal(x)) as(x, "diagonalMatrix") else x      if(isS4(x) && isDiagonal(x)) as(x, "diagonalMatrix") else x
# Line 345  Line 348 
348                           j = "index", drop = "logical"), subDiag)                           j = "index", drop = "logical"), subDiag)
349  setMethod("[", signature(x = "diagonalMatrix", i = "index",  setMethod("[", signature(x = "diagonalMatrix", i = "index",
350                          j = "missing", drop = "logical"),                          j = "missing", drop = "logical"),
351            function(x, i, j, ..., drop) subDiag(x, i=i, drop=drop))            function(x, i, j, ..., drop) {
352                  na <- nargs()
353                  Matrix.msg("diag[i,m,l] : nargs()=", na, .M.level = 2)
354                  if(na == 4)
355                       subDiag(x, i=i, , drop=drop)
356                  else subDiag(x, i=i,   drop=drop)
357              })
358  setMethod("[", signature(x = "diagonalMatrix", i = "missing",  setMethod("[", signature(x = "diagonalMatrix", i = "missing",
359                           j = "index", drop = "logical"),                           j = "index", drop = "logical"),
360            function(x, i, j, ..., drop) subDiag(x, j=j, drop=drop))            function(x, i, j, ..., drop) subDiag(x, j=j, drop=drop))
# Line 355  Line 364 
364  ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch  ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch
365  ## Only(?) current bug:  x[i] <- value  is wrong when  i is *vector*  ## Only(?) current bug:  x[i] <- value  is wrong when  i is *vector*
366  replDiag <- function(x, i, j, ..., value) {  replDiag <- function(x, i, j, ..., value) {
367      x <- as(x, "sparseMatrix")      x <- as(x, "TsparseMatrix")
368      if(missing(i))      if(missing(i))
369          x[, j] <- value          x[, j] <- value
370      else if(missing(j)) { ##  x[i , ] <- v  *OR*   x[i] <- v      else if(missing(j)) { ##  x[i , ] <- v  *OR*   x[i] <- v
# Line 400  Line 409 
409                               x@x[ii] <- value                               x@x[ii] <- value
410                               x                               x
411                           } else { ## no longer diagonal, but remain sparse:                           } else { ## no longer diagonal, but remain sparse:
412                               x <- as(x, "sparseMatrix")                               x <- as(x, "TsparseMatrix")
413                               x[i] <- value                               x[i] <- value
414                               x                               x
415                           }                           }
# Line 417  Line 426 
426                   function(x,i,j, ..., value) replDiag(x, j=j, value=value))                   function(x,i,j, ..., value) replDiag(x, j=j, value=value))
427    
428  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",
429                                  value = "scarceMatrix"),                                  value = "sparseMatrix"),
430                   function (x, i, j, ..., value)                   function (x, i, j, ..., value)
431                   callGeneric(x=x, , j=j, value = as(value, "sparseVector")))                   callGeneric(x=x, , j=j, value = as(value, "sparseVector")))
432  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "missing",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "missing",
433                                  value = "scarceMatrix"),                                  value = "sparseMatrix"),
434                   function (x, i, j, ..., value)                   function (x, i, j, ..., value)
435                   callGeneric(x=x, i=i, , value = as(value, "sparseVector")))                   callGeneric(x=x, i=i, , value = as(value, "sparseVector")))
436  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "index",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "index",
437                                  value = "scarceMatrix"),                                  value = "sparseMatrix"),
438                   function (x, i, j, ..., value)                   function (x, i, j, ..., value)
439                   callGeneric(x=x, i=i, j=j, value = as(value, "sparseVector")))                   callGeneric(x=x, i=i, j=j, value = as(value, "sparseVector")))
440    
# Line 591  Line 600 
600  ##           })  ##           })
601    
602  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
603            function(x, y = NULL) crossprod(as(x, "sparseMatrix"), y))            function(x, y = NULL) crossprod(as(x, "TsparseMatrix"), y))
604    
605  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
606            function(x, y = NULL) crossprod(x, as(y, "sparseMatrix")))            function(x, y = NULL) crossprod(x, as(y, "TsparseMatrix")))
607    
608  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
609            function(x, y = NULL) tcrossprod(as(x, "sparseMatrix"), y))            function(x, y = NULL) tcrossprod(as(x, "TsparseMatrix"), y))
610    
611  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
612            function(x, y = NULL) tcrossprod(x, as(y, "sparseMatrix")))            function(x, y = NULL) tcrossprod(x, as(y, "TsparseMatrix")))
613    
614    
615  ## FIXME?: In theory, this can be done *FASTER*, in some cases, via tapply1()  ## FIXME?: In theory, this can be done *FASTER*, in some cases, via tapply1()
616  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
617            function(x, y) as(x, "sparseMatrix") %*% y)            function(x, y) as(x, "TsparseMatrix") %*% y)
618  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
619            function(x, y) x %*% as(y, "sparseMatrix"))            function(x, y) x %*% as(y, "TsparseMatrix"))
620  ## NB: The previous is *not* triggering for  "ddi" o "dgC" (= distance 3)  ## NB: The previous is *not* triggering for  "ddi" o "dgC" (= distance 3)
621  ##     since there's a "ddense" o "Csparse" at dist. 2 => triggers first.  ##     since there's a "ddense" o "Csparse" at dist. 2 => triggers first.
622  ## ==> do this:  ## ==> do this:
# Line 673  Line 682 
682          isNum <- (is.numeric(r) || is.numeric(r00))          isNum <- (is.numeric(r) || is.numeric(r00))
683          isLog <- (is.logical(r) || is.logical(r00))          isLog <- (is.logical(r) || is.logical(r00))
684    
685          if(getOption("verbose"))          Matrix.msg("exploding   <diag>  o  <diag>  into dense matrix")
             message("exploding  <diag>  o  <diag>  into dense matrix")  
686          d <- e1@Dim          d <- e1@Dim
687          n <- d[1]          n <- d[1]
688          stopifnot(length(r) == n)          stopifnot(length(r) == n)
# Line 897  Line 905 
905      setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),      setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),
906                function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "l")))                function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "l")))
907  }  }
908  ## This should *not* dispatch to <dense> methods (in ./Ops.R ), as  
909  ##  FALSE & <anything> |-> FALSE : hence result should be diagonal:  ## Direct subclasses of "denseMatrix": currently ddenseMatrix, ldense... :
910  for(cl in diCls) {  dense.subCl <- local({ dM.scl <- getClass("denseMatrix")@subclasses
911      setMethod("&", signature(e1 = cl, e2 = "ANY"),                         names(dM.scl)[sapply(dM.scl, slot, "distance") == 1] })
912                function(e1,e2) e1 & as(e2,"Matrix"))  for(DI in diCls) {
913      setMethod("&", signature(e1 = "ANY", e2 = cl),      for(c2 in c(dense.subCl, "Matrix")) {
914                function(e1,e2) as(e1,"Matrix") & e2)          for(Fun in c("*", "^", "&")) {
915      for(c2 in c("denseMatrix", "Matrix")) {              setMethod(Fun, signature(e1 = DI, e2 = c2),
916          setMethod("&", signature(e1 = cl, e2 = c2),                        function(e1,e2) callGeneric(e1, Diagonal(x = diag(e2))))
917                    function(e1,e2) e1 & Diagonal(x = diag(e2)))              setMethod(Fun, signature(e1 = c2, e2 = DI),
918          setMethod("&", signature(e1 = c2, e2 = cl),                        function(e1,e2) callGeneric(Diagonal(x = diag(e1)), e2))
919                    function(e1,e2) Diagonal(x = diag(e1)) & e2)          }
920            ## NB: This arguably implicitly uses  0/0 :== 0  to keep diagonality
921            for(Fun in c("%%", "%/%", "/")) {
922                setMethod(Fun, signature(e1 = DI, e2 = c2),
923                          function(e1,e2) callGeneric(e1, Diagonal(x = diag(e2))))
924            }
925      }      }
926  }  }
927    
928    
929  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"
930  ### ----------  any, all: separately here  ### ----------   the last 4: separately here
931  for(cl in diCls) {  for(cl in diCls) {
932  setMethod("any", cl,  setMethod("any", cl,
933            function (x, ..., na.rm) {            function (x, ..., na.rm) {
934                if(any(x@Dim == 0)) FALSE                if(any(x@Dim == 0)) FALSE
935                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)
936            })            })
937  setMethod("all",  cl, function (x, ..., na.rm) any(x@Dim == 0))  setMethod("all",  cl, function (x, ..., na.rm) {
938  setMethod("prod", cl, function (x, ..., na.rm) as.numeric(any(x@Dim == 0)))      n <- x@Dim[1]
939        if(n >= 2) FALSE
940        else if(n == 0 || x@diag == "U") TRUE
941        else all(x@x, ..., na.rm = na.rm)
942    })
943    setMethod("prod", cl, function (x, ..., na.rm) {
944        n <- x@Dim[1]
945        if(n >= 2) 0
946        else if(n == 0 || x@diag == "U") 1
947        else ## n == 1, diag = "N" :
948            prod(x@x, ..., na.rm = na.rm)
949    })
950    
951  setMethod("sum", cl,  setMethod("sum", cl,
952            function(x, ..., na.rm) {            function(x, ..., na.rm) {

Legend:
Removed from v.2237  
changed lines
  Added in v.2363

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