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 2417, Fri Jul 10 16:17:28 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 586  Line 595 
595  ##        function(x, y = NULL) {  ##        function(x, y = NULL) {
596  ##           })  ##           })
597    
598  ## setMethod("tcrossprod", signature(x = "denseMatrix", y = "diagonalMatrix"),  Cspdiagprod <- function(x, y) {
599  ##        function(x, y = NULL) {      dx <- dim(x)
600  ##           })      dy <- dim(y)
601        if(dx[2] != dy[1]) stop("non-matching dimensions")
602        ind <- rep.int(seq_len(dx[2]), x@p[-1] - x@p[-dx[2]-1L])
603        if(y@diag == "N")
604            x@x <- x@x * y@x[ind]
605        x
606    }
607    
608    diagCspprod <- function(x, y) {
609        dx <- dim(x)
610        dy <- dim(y)
611        if(dx[2] != dy[1]) stop("non-matching dimensions")
612        if(x@diag == "N")
613            y@x <- y@x * x@x[y@i + 1L]
614        y
615    }
616    
617    setMethod("crossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
618              function(x, y = NULL) diagCspprod(x, y))
619    
620  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
621            function(x, y = NULL) crossprod(as(x, "sparseMatrix"), y))            function(x, y = NULL) diagCspprod(x, as(y, "CsparseMatrix")))
622    
623    ## Prefer calling diagCspprod to Cspdiagprod if going to transpose anyway
624    ##  x'y == (y'x)'
625    setMethod("crossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
626              function(x, y = NULL) t(diagCspprod(y, x)))
627    
628  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
629            function(x, y = NULL) crossprod(x, as(y, "sparseMatrix")))            function(x, y = NULL) t(diagCspprod(y, as(x, "Csparsematrix"))))
630    
631    setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
632              function(x, y = NULL) diagCspprod(x, t(y)))
633    
634  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
635            function(x, y = NULL) tcrossprod(as(x, "sparseMatrix"), y))            function(x, y = NULL) diagCspprod(x, t(as(y, "CsparseMatrix"))))
636    
637    setMethod("tcrossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
638              function(x, y = NULL) Cspdiagprod(x, y))
639    
640  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
641            function(x, y = NULL) tcrossprod(x, as(y, "sparseMatrix")))            function(x, y = NULL) Cspdiagprod(as(x, "CsparseMatrix"), y))
642    
643    setMethod("%*%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
644              function(x, y) diagCspprod(x, y))
645    
 ## FIXME?: In theory, this can be done *FASTER*, in some cases, via tapply1()  
646  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
647            function(x, y) as(x, "sparseMatrix") %*% y)            function(x, y) diagCspprod(as(x, "CsparseMatrix"), y))
648    
649  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
650            function(x, y) x %*% as(y, "sparseMatrix"))            function(x, y) Cspdiagprod(as(x, "CsparseMatrix"), y))
651  ## 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)  
652  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
653            function(x, y) x %*% as(y, "CsparseMatrix"))            function(x, y) Cspdiagprod(x, y))
654  ## NB: this is *not* needed for Tsparse & Rsparse  
655  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*
656  ##       do indeed work by going through sparse (and *not* ddense)!  ##       do indeed work by going through sparse (and *not* ddense)!
657    
   
   
658  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),
659            function(a, b, ...) {            function(a, b, ...) {
660                a@x <- 1/ a@x                a@x <- 1/ a@x
# Line 673  Line 707 
707          isNum <- (is.numeric(r) || is.numeric(r00))          isNum <- (is.numeric(r) || is.numeric(r00))
708          isLog <- (is.logical(r) || is.logical(r00))          isLog <- (is.logical(r) || is.logical(r00))
709    
710          if(getOption("verbose"))          Matrix.msg("exploding   <diag>  o  <diag>  into dense matrix")
             message("exploding  <diag>  o  <diag>  into dense matrix")  
711          d <- e1@Dim          d <- e1@Dim
712          n <- d[1]          n <- d[1]
713          stopifnot(length(r) == n)          stopifnot(length(r) == n)
# Line 897  Line 930 
930      setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),      setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),
931                function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "l")))                function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "l")))
932  }  }
933  ## This should *not* dispatch to <dense> methods (in ./Ops.R ), as  
934  ##  FALSE & <anything> |-> FALSE : hence result should be diagonal:  ## Direct subclasses of "denseMatrix": currently ddenseMatrix, ldense... :
935  for(cl in diCls) {  dense.subCl <- local({ dM.scl <- getClass("denseMatrix")@subclasses
936      setMethod("&", signature(e1 = cl, e2 = "ANY"),                         names(dM.scl)[sapply(dM.scl, slot, "distance") == 1] })
937                function(e1,e2) e1 & as(e2,"Matrix"))  for(DI in diCls) {
938      setMethod("&", signature(e1 = "ANY", e2 = cl),      for(c2 in c(dense.subCl, "Matrix")) {
939                function(e1,e2) as(e1,"Matrix") & e2)          for(Fun in c("*", "^", "&")) {
940      for(c2 in c("denseMatrix", "Matrix")) {              setMethod(Fun, signature(e1 = DI, e2 = c2),
941          setMethod("&", signature(e1 = cl, e2 = c2),                        function(e1,e2) callGeneric(e1, Diagonal(x = diag(e2))))
942                    function(e1,e2) e1 & Diagonal(x = diag(e2)))              setMethod(Fun, signature(e1 = c2, e2 = DI),
943          setMethod("&", signature(e1 = c2, e2 = cl),                        function(e1,e2) callGeneric(Diagonal(x = diag(e1)), e2))
944                    function(e1,e2) Diagonal(x = diag(e1)) & e2)          }
945            ## NB: This arguably implicitly uses  0/0 :== 0  to keep diagonality
946            for(Fun in c("%%", "%/%", "/")) {
947                setMethod(Fun, signature(e1 = DI, e2 = c2),
948                          function(e1,e2) callGeneric(e1, Diagonal(x = diag(e2))))
949            }
950      }      }
951  }  }
952    
953    
954  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"
955  ### ----------  any, all: separately here  ### ----------   the last 4: separately here
956  for(cl in diCls) {  for(cl in diCls) {
957  setMethod("any", cl,  setMethod("any", cl,
958            function (x, ..., na.rm) {            function (x, ..., na.rm) {
959                if(any(x@Dim == 0)) FALSE                if(any(x@Dim == 0)) FALSE
960                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)
961            })            })
962  setMethod("all",  cl, function (x, ..., na.rm) any(x@Dim == 0))  setMethod("all",  cl, function (x, ..., na.rm) {
963  setMethod("prod", cl, function (x, ..., na.rm) as.numeric(any(x@Dim == 0)))      n <- x@Dim[1]
964        if(n >= 2) FALSE
965        else if(n == 0 || x@diag == "U") TRUE
966        else all(x@x, ..., na.rm = na.rm)
967    })
968    setMethod("prod", cl, function (x, ..., na.rm) {
969        n <- x@Dim[1]
970        if(n >= 2) 0
971        else if(n == 0 || x@diag == "U") 1
972        else ## n == 1, diag = "N" :
973            prod(x@x, ..., na.rm = na.rm)
974    })
975    
976  setMethod("sum", cl,  setMethod("sum", cl,
977            function(x, ..., na.rm) {            function(x, ..., na.rm) {

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

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