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 2341, Mon Mar 2 17:53:13 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))
# Line 355  Line 358 
358  ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch  ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch
359  ## Only(?) current bug:  x[i] <- value  is wrong when  i is *vector*  ## Only(?) current bug:  x[i] <- value  is wrong when  i is *vector*
360  replDiag <- function(x, i, j, ..., value) {  replDiag <- function(x, i, j, ..., value) {
361      x <- as(x, "sparseMatrix")      x <- as(x, "TsparseMatrix")
362      if(missing(i))      if(missing(i))
363          x[, j] <- value          x[, j] <- value
364      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 403 
403                               x@x[ii] <- value                               x@x[ii] <- value
404                               x                               x
405                           } else { ## no longer diagonal, but remain sparse:                           } else { ## no longer diagonal, but remain sparse:
406                               x <- as(x, "sparseMatrix")                               x <- as(x, "TsparseMatrix")
407                               x[i] <- value                               x[i] <- value
408                               x                               x
409                           }                           }
# Line 417  Line 420 
420                   function(x,i,j, ..., value) replDiag(x, j=j, value=value))                   function(x,i,j, ..., value) replDiag(x, j=j, value=value))
421    
422  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",
423                                  value = "scarceMatrix"),                                  value = "sparseMatrix"),
424                   function (x, i, j, ..., value)                   function (x, i, j, ..., value)
425                   callGeneric(x=x, , j=j, value = as(value, "sparseVector")))                   callGeneric(x=x, , j=j, value = as(value, "sparseVector")))
426  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "missing",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "missing",
427                                  value = "scarceMatrix"),                                  value = "sparseMatrix"),
428                   function (x, i, j, ..., value)                   function (x, i, j, ..., value)
429                   callGeneric(x=x, i=i, , value = as(value, "sparseVector")))                   callGeneric(x=x, i=i, , value = as(value, "sparseVector")))
430  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "index",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "index",
431                                  value = "scarceMatrix"),                                  value = "sparseMatrix"),
432                   function (x, i, j, ..., value)                   function (x, i, j, ..., value)
433                   callGeneric(x=x, i=i, j=j, value = as(value, "sparseVector")))                   callGeneric(x=x, i=i, j=j, value = as(value, "sparseVector")))
434    
# Line 591  Line 594 
594  ##           })  ##           })
595    
596  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
597            function(x, y = NULL) crossprod(as(x, "sparseMatrix"), y))            function(x, y = NULL) crossprod(as(x, "TsparseMatrix"), y))
598    
599  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
600            function(x, y = NULL) crossprod(x, as(y, "sparseMatrix")))            function(x, y = NULL) crossprod(x, as(y, "TsparseMatrix")))
601    
602  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
603            function(x, y = NULL) tcrossprod(as(x, "sparseMatrix"), y))            function(x, y = NULL) tcrossprod(as(x, "TsparseMatrix"), y))
604    
605  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
606            function(x, y = NULL) tcrossprod(x, as(y, "sparseMatrix")))            function(x, y = NULL) tcrossprod(x, as(y, "TsparseMatrix")))
607    
608    
609  ## 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()
610  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
611            function(x, y) as(x, "sparseMatrix") %*% y)            function(x, y) as(x, "TsparseMatrix") %*% y)
612  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
613            function(x, y) x %*% as(y, "sparseMatrix"))            function(x, y) x %*% as(y, "TsparseMatrix"))
614  ## NB: The previous is *not* triggering for  "ddi" o "dgC" (= distance 3)  ## NB: The previous is *not* triggering for  "ddi" o "dgC" (= distance 3)
615  ##     since there's a "ddense" o "Csparse" at dist. 2 => triggers first.  ##     since there's a "ddense" o "Csparse" at dist. 2 => triggers first.
616  ## ==> do this:  ## ==> do this:
# Line 897  Line 900 
900      setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),      setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),
901                function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "l")))                function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "l")))
902  }  }
903  ## This should *not* dispatch to <dense> methods (in ./Ops.R ), as  
904  ##  FALSE & <anything> |-> FALSE : hence result should be diagonal:  ## Direct subclasses of "denseMatrix": currently ddenseMatrix, ldense... :
905  for(cl in diCls) {  dense.subCl <- local({ dM.scl <- getClass("denseMatrix")@subclasses
906      setMethod("&", signature(e1 = cl, e2 = "ANY"),                         names(dM.scl)[sapply(dM.scl, slot, "distance") == 1] })
907                function(e1,e2) e1 & as(e2,"Matrix"))  for(DI in diCls) {
908      setMethod("&", signature(e1 = "ANY", e2 = cl),      for(c2 in c(dense.subCl, "Matrix")) {
909                function(e1,e2) as(e1,"Matrix") & e2)          for(Fun in c("*", "^", "&")) {
910      for(c2 in c("denseMatrix", "Matrix")) {              setMethod(Fun, signature(e1 = DI, e2 = c2),
911          setMethod("&", signature(e1 = cl, e2 = c2),                        function(e1,e2) callGeneric(e1, Diagonal(x = diag(e2))))
912                    function(e1,e2) e1 & Diagonal(x = diag(e2)))              setMethod(Fun, signature(e1 = c2, e2 = DI),
913          setMethod("&", signature(e1 = c2, e2 = cl),                        function(e1,e2) callGeneric(Diagonal(x = diag(e1)), e2))
914                    function(e1,e2) Diagonal(x = diag(e1)) & e2)          }
915            ## NB: This arguably implicitly uses  0/0 :== 0  to keep diagonality
916            for(Fun in c("%%", "%/%", "/")) {
917                setMethod(Fun, signature(e1 = DI, e2 = c2),
918                          function(e1,e2) callGeneric(e1, Diagonal(x = diag(e2))))
919            }
920      }      }
921  }  }
922    
# Line 921  Line 929 
929                if(any(x@Dim == 0)) FALSE                if(any(x@Dim == 0)) FALSE
930                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)
931            })            })
932  setMethod("all",  cl, function (x, ..., na.rm) any(x@Dim == 0))  setMethod("all",  cl, function (x, ..., na.rm) {
933  setMethod("prod", cl, function (x, ..., na.rm) as.numeric(any(x@Dim == 0)))      n <- x@Dim[1]
934        if(n >= 2) FALSE
935        else if(n == 0 || x@diag == "U") TRUE
936        else all(x@x, ..., na.rm = na.rm)
937    })
938    setMethod("prod", cl, function (x, ..., na.rm) {
939        n <- x@Dim[1]
940        if(n >= 2) 0
941        else if(n == 0 || x@diag == "U") 1
942        else ## n == 1, diag = "N" :
943            prod(x@x, ..., na.rm = na.rm)
944    })
945    
946  setMethod("sum", cl,  setMethod("sum", cl,
947            function(x, ..., na.rm) {            function(x, ..., na.rm) {

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

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