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 2496, Sat Nov 14 17:24:42 2009 UTC
# Line 45  Line 45 
45              storage.mode(x) <- "double"              storage.mode(x) <- "double"
46              "d"              "d"
47          }          }
48      new(paste(kind, shape, "CMatrix", sep=''),      ii <- if(n) 0:(n - 1L) else integer(0)
49          Dim = c(n,n), x = x, uplo = uplo,      if(shape == "g")
50          i = if(n) 0:(n - 1L) else integer(0), p = 0:n)          new(paste0(kind, "gCMatrix"), Dim = c(n,n),
51                x = x, i = ii, p = 0:n)
52        else new(paste0(kind, shape, "CMatrix"), Dim = c(n,n), uplo = uplo,
53                 x = x, i = ii, p = 0:n)
54  }  }
55    
56  ## Pkg 'spdep' had (relatively slow) versions of this as_dsCMatrix_I()  ## Pkg 'spdep' had (relatively slow) versions of this as_dsCMatrix_I()
# Line 214  Line 217 
217  ## ddi*:  ## ddi*:
218  diag2tT <- function(from) .diag2tT(from, "U", "d")  diag2tT <- function(from) .diag2tT(from, "U", "d")
219  setAs("ddiMatrix", "triangularMatrix", diag2tT)  setAs("ddiMatrix", "triangularMatrix", diag2tT)
220  setAs("ddiMatrix", "sparseMatrix", diag2tT)  ##_no_longer_ setAs("ddiMatrix", "sparseMatrix", diag2tT)
221  ## needed too (otherwise <dense> -> Tsparse is taken):  ## needed too (otherwise <dense> -> Tsparse is taken):
222  setAs("ddiMatrix", "TsparseMatrix", diag2tT)  setAs("ddiMatrix", "TsparseMatrix", diag2tT)
223    setAs("ddiMatrix", "dsparseMatrix", diag2tT)
224  setAs("ddiMatrix", "CsparseMatrix",  setAs("ddiMatrix", "CsparseMatrix",
225        function(from) as(.diag2tT(from, "U", "d"), "CsparseMatrix"))        function(from) as(.diag2tT(from, "U", "d"), "CsparseMatrix"))
226  setAs("ddiMatrix", "symmetricMatrix",  setAs("ddiMatrix", "symmetricMatrix",
# Line 225  Line 229 
229  ## ldi*:  ## ldi*:
230  diag2tT <- function(from) .diag2tT(from, "U", "l")  diag2tT <- function(from) .diag2tT(from, "U", "l")
231  setAs("ldiMatrix", "triangularMatrix", diag2tT)  setAs("ldiMatrix", "triangularMatrix", diag2tT)
232  setAs("ldiMatrix", "sparseMatrix", diag2tT)  ##_no_longer_ setAs("ldiMatrix", "sparseMatrix", diag2tT)
233  ## needed too (otherwise <dense> -> Tsparse is taken):  ## needed too (otherwise <dense> -> Tsparse is taken):
234  setAs("ldiMatrix", "TsparseMatrix", diag2tT)  setAs("ldiMatrix", "TsparseMatrix", diag2tT)
235    setAs("ldiMatrix", "lsparseMatrix", diag2tT)
236  setAs("ldiMatrix", "CsparseMatrix",  setAs("ldiMatrix", "CsparseMatrix",
237        function(from) as(.diag2tT(from, "U", "l"), "CsparseMatrix"))        function(from) as(.diag2tT(from, "U", "l"), "CsparseMatrix"))
238  setAs("ldiMatrix", "symmetricMatrix",  setAs("ldiMatrix", "symmetricMatrix",
# Line 272  Line 277 
277  setAs("diagonalMatrix", "generalMatrix", # prefer sparse:  setAs("diagonalMatrix", "generalMatrix", # prefer sparse:
278        function(from) as(as(from, "CsparseMatrix"), "generalMatrix"))        function(from) as(as(from, "CsparseMatrix"), "generalMatrix"))
279    
280    setAs("diagonalMatrix", "denseMatrix",
281          function(from) as(as(from, "CsparseMatrix"), "denseMatrix"))
282    
283  .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
284    
285  .diag.2N <- function(m) {  .diag.2N <- function(m) {
# Line 292  Line 300 
300            d <- dim(from)            d <- dim(from)
301            if(d[1] != (n <- d[2])) stop("non-square matrix")            if(d[1] != (n <- d[2])) stop("non-square matrix")
302            if(any(from[row(from) != col(from)] != 0))            if(any(from[row(from) != col(from)] != 0))
303                stop("has non-zero off-diagonal entries")                stop("matrix with non-zero off-diagonals cannot be coerced to diagonalMatrix")
304            x <- diag(from)            x <- diag(from)
305            if(is.logical(x)) {            if(is.logical(x)) {
306                cl <- "ldiMatrix"                cl <- "ldiMatrix"
307                uni <- all(x)                uni <- allTrue(x) ## uni := {is it unit-diagonal ?}
308            } else {            } else {
309                cl <- "ddiMatrix"                cl <- "ddiMatrix"
310                uni <- all(x == 1)                uni <- allTrue(x == 1)
311                storage.mode(x) <- "double"                storage.mode(x) <- "double"
312            } ## TODO: complex            } ## TODO: complex
313            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 324 
324            x <- diag(from)            x <- diag(from)
325            if(is.logical(x)) {            if(is.logical(x)) {
326                cl <- "ldiMatrix"                cl <- "ldiMatrix"
327                uni <- all(x)                uni <- allTrue(x)
328            } else {            } else {
329                cl <- "ddiMatrix"                cl <- "ddiMatrix"
330                uni <- all(x == 1)                uni <- allTrue(x == 1)
331                storage.mode(x) <- "double"                storage.mode(x) <- "double"
332            }            } ## TODO: complex
333            new(cl, Dim = c(n,n), diag = if(uni) "U" else "N",            new(cl, Dim = c(n,n), diag = if(uni) "U" else "N",
334                x = if(uni) x[FALSE] else x)                x = if(uni) x[FALSE] else x)
335        })        })
# Line 331  Line 339 
339            function(x = 1, nrow, ncol) .diag.x(x))            function(x = 1, nrow, ncol) .diag.x(x))
340    
341  subDiag <- function(x, i, j, ..., drop) {  subDiag <- function(x, i, j, ..., drop) {
342      x <- as(x, "sparseMatrix")      x <- as(x, "TsparseMatrix")
343      x <- if(missing(i))      x <- if(missing(i))
344          x[, j, drop=drop]          x[, j, drop=drop]
345      else if(missing(j))      else if(missing(j))
346          x[i, , drop=drop]          if(nargs() == 4) x[i, , drop=drop] else x[i, drop=drop]
347      else      else
348          x[i,j, drop=drop]          x[i,j, drop=drop]
349      if(isS4(x) && isDiagonal(x)) as(x, "diagonalMatrix") else x      if(isS4(x) && isDiagonal(x)) as(x, "diagonalMatrix") else x
# Line 345  Line 353 
353                           j = "index", drop = "logical"), subDiag)                           j = "index", drop = "logical"), subDiag)
354  setMethod("[", signature(x = "diagonalMatrix", i = "index",  setMethod("[", signature(x = "diagonalMatrix", i = "index",
355                          j = "missing", drop = "logical"),                          j = "missing", drop = "logical"),
356            function(x, i, j, ..., drop) subDiag(x, i=i, drop=drop))            function(x, i, j, ..., drop) {
357                  na <- nargs()
358                  Matrix.msg("diag[i,m,l] : nargs()=", na, .M.level = 2)
359                  if(na == 4)
360                       subDiag(x, i=i, , drop=drop)
361                  else subDiag(x, i=i,   drop=drop)
362              })
363  setMethod("[", signature(x = "diagonalMatrix", i = "missing",  setMethod("[", signature(x = "diagonalMatrix", i = "missing",
364                           j = "index", drop = "logical"),                           j = "index", drop = "logical"),
365            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 369 
369  ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch  ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch
370  ## Only(?) current bug:  x[i] <- value  is wrong when  i is *vector*  ## Only(?) current bug:  x[i] <- value  is wrong when  i is *vector*
371  replDiag <- function(x, i, j, ..., value) {  replDiag <- function(x, i, j, ..., value) {
372      x <- as(x, "sparseMatrix")      x <- as(x, "TsparseMatrix")
373      if(missing(i))      if(missing(i))
374          x[, j] <- value          x[, j] <- value
375      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 414 
414                               x@x[ii] <- value                               x@x[ii] <- value
415                               x                               x
416                           } else { ## no longer diagonal, but remain sparse:                           } else { ## no longer diagonal, but remain sparse:
417                               x <- as(x, "sparseMatrix")                               x <- as(x, "TsparseMatrix")
418                               x[i] <- value                               x[i] <- value
419                               x                               x
420                           }                           }
# Line 417  Line 431 
431                   function(x,i,j, ..., value) replDiag(x, j=j, value=value))                   function(x,i,j, ..., value) replDiag(x, j=j, value=value))
432    
433  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",
434                                  value = "scarceMatrix"),                                  value = "sparseMatrix"),
435                   function (x, i, j, ..., value)                   function (x, i, j, ..., value)
436                   callGeneric(x=x, , j=j, value = as(value, "sparseVector")))                   callGeneric(x=x, , j=j, value = as(value, "sparseVector")))
437  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "missing",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "missing",
438                                  value = "scarceMatrix"),                                  value = "sparseMatrix"),
439                   function (x, i, j, ..., value)                   function (x, i, j, ..., value)
440                   callGeneric(x=x, i=i, , value = as(value, "sparseVector")))                   callGeneric(x=x, i=i, , value = as(value, "sparseVector")))
441  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "index",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "index",
442                                  value = "scarceMatrix"),                                  value = "sparseMatrix"),
443                   function (x, i, j, ..., value)                   function (x, i, j, ..., value)
444                   callGeneric(x=x, i=i, j=j, value = as(value, "sparseVector")))                   callGeneric(x=x, i=i, j=j, value = as(value, "sparseVector")))
445    
# Line 586  Line 600 
600  ##        function(x, y = NULL) {  ##        function(x, y = NULL) {
601  ##           })  ##           })
602    
603  ## setMethod("tcrossprod", signature(x = "denseMatrix", y = "diagonalMatrix"),  Cspdiagprod <- function(x, y) {
604  ##        function(x, y = NULL) {      dx <- dim(x <- .Call(Csparse_diagU2N, x))
605  ##           })      dy <- dim(y)
606        if(dx[2] != dy[1]) stop("non-matching dimensions")
607        ind <- rep.int(seq_len(dx[2]), x@p[-1] - x@p[-dx[2]-1L])
608        if(y@diag == "N")
609            x@x <- x@x * y@x[ind]
610        x
611    }
612    
613    diagCspprod <- function(x, y) {
614        dx <- dim(x)
615        dy <- dim(y <- .Call(Csparse_diagU2N, y))
616        if(dx[2] != dy[1]) stop("non-matching dimensions")
617        if(x@diag == "N")
618            y@x <- y@x * x@x[y@i + 1L]
619        y
620    }
621    
622    setMethod("crossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
623              function(x, y = NULL) diagCspprod(x, y))
624    
625  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
626            function(x, y = NULL) crossprod(as(x, "sparseMatrix"), y))            function(x, y = NULL) diagCspprod(x, as(y, "CsparseMatrix")))
627    
628    ## Prefer calling diagCspprod to Cspdiagprod if going to transpose anyway
629    ##  x'y == (y'x)'
630    setMethod("crossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
631              function(x, y = NULL) t(diagCspprod(y, x)))
632    
633  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
634            function(x, y = NULL) crossprod(x, as(y, "sparseMatrix")))            function(x, y = NULL) t(diagCspprod(y, as(x, "Csparsematrix"))))
635    
636    setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
637              function(x, y = NULL) diagCspprod(x, t(y)))
638    
639  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
640            function(x, y = NULL) tcrossprod(as(x, "sparseMatrix"), y))            function(x, y = NULL) diagCspprod(x, t(as(y, "CsparseMatrix"))))
641    
642    setMethod("tcrossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
643              function(x, y = NULL) Cspdiagprod(x, y))
644    
645  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
646            function(x, y = NULL) tcrossprod(x, as(y, "sparseMatrix")))            function(x, y = NULL) Cspdiagprod(as(x, "CsparseMatrix"), y))
647    
648    setMethod("%*%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
649              function(x, y) diagCspprod(x, y))
650    
 ## FIXME?: In theory, this can be done *FASTER*, in some cases, via tapply1()  
651  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
652            function(x, y) as(x, "sparseMatrix") %*% y)            function(x, y) diagCspprod(as(x, "CsparseMatrix"), y))
653    
654  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
655            function(x, y) x %*% as(y, "sparseMatrix"))            function(x, y) Cspdiagprod(as(x, "CsparseMatrix"), y))
656  ## 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)  
657  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
658            function(x, y) x %*% as(y, "CsparseMatrix"))            function(x, y) Cspdiagprod(x, y))
659  ## NB: this is *not* needed for Tsparse & Rsparse  
660  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*
661  ##       do indeed work by going through sparse (and *not* ddense)!  ##       do indeed work by going through sparse (and *not* ddense)!
662    
   
   
663  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),
664            function(a, b, ...) {            function(a, b, ...) {
665                a@x <- 1/ a@x                a@x <- 1/ a@x
# Line 670  Line 709 
709          .diag.2N(e1)          .diag.2N(e1)
710      }      }
711      else { ## result not diagonal, but at least symmetric:      else { ## result not diagonal, but at least symmetric:
712            ## e.g., m == m
713          isNum <- (is.numeric(r) || is.numeric(r00))          isNum <- (is.numeric(r) || is.numeric(r00))
714          isLog <- (is.logical(r) || is.logical(r00))          isLog <- (is.logical(r) || is.logical(r00))
715            Matrix.msg("exploding <diag> o <diag> into dense matrix", .M.level = 2)
         if(getOption("verbose"))  
             message("exploding  <diag>  o  <diag>  into dense matrix")  
716          d <- e1@Dim          d <- e1@Dim
717          n <- d[1]          n <- d[1]
718          stopifnot(length(r) == n)          stopifnot(length(r) == n)
# Line 702  Line 740 
740              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)
741  }  }
742    
743  ## FIXME:    diagonal  o  triangular  |-->  triangular  ## diagonal  o  triangular  |-->  triangular
744  ## -----     diagonal  o  symmetric   |-->  symmetric  ## diagonal  o  symmetric   |-->  symmetric
745  ##    {also when other is sparse: do these "here" --  ##    {also when other is sparse: do these "here" --
746  ##     before conversion to sparse, since that loses "diagonality"}  ##     before conversion to sparse, since that loses "diagonality"}
747    
# Line 888  Line 926 
926  for(other in c("ANY", "Matrix", "dMatrix")) {  for(other in c("ANY", "Matrix", "dMatrix")) {
927      ## ddi*:      ## ddi*:
928      setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other),      setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other),
929                function(e1,e2) callGeneric(diag2tT.u(e1,e2, "d"), e2))                function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2))
930      setMethod("Ops", signature(e1 = other, e2 = "ddiMatrix"),      setMethod("Ops", signature(e1 = other, e2 = "ddiMatrix"),
931                function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "d")))                function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "d")))
932      ## ldi*:      ## ldi*:
933      setMethod("Ops", signature(e1 = "ldiMatrix", e2 = other),      setMethod("Ops", signature(e1 = "ldiMatrix", e2 = other),
934                function(e1,e2) callGeneric(diag2tT.u(e1,e2, "l"), e2))                function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2))
935      setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),      setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),
936                function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "l")))                function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "l")))
937    }
938    
939    ## Direct subclasses of "denseMatrix": currently ddenseMatrix, ldense... :
940    dense.subCl <- local({ dM.scl <- getClass("denseMatrix")@subclasses
941                           names(dM.scl)[sapply(dM.scl, slot, "distance") == 1] })
942    for(DI in diCls) {
943        for(c2 in c(dense.subCl, "Matrix")) {
944            for(Fun in c("*", "^", "&")) {
945                setMethod(Fun, signature(e1 = DI, e2 = c2),
946                          function(e1,e2) callGeneric(e1, Diagonal(x = diag(e2))))
947                setMethod(Fun, signature(e1 = c2, e2 = DI),
948                          function(e1,e2) callGeneric(Diagonal(x = diag(e1)), e2))
949            }
950            ## NB: This arguably implicitly uses  0/0 :== 0  to keep diagonality
951            for(Fun in c("%%", "%/%", "/")) {
952                setMethod(Fun, signature(e1 = DI, e2 = c2),
953                          function(e1,e2) callGeneric(e1, Diagonal(x = diag(e2))))
954  }  }
 ## This should *not* dispatch to <dense> methods (in ./Ops.R ), as  
 ##  FALSE & <anything> |-> FALSE : hence result should be diagonal:  
 for(cl in diCls) {  
     setMethod("&", signature(e1 = cl, e2 = "ANY"),  
               function(e1,e2) e1 & as(e2,"Matrix"))  
     setMethod("&", signature(e1 = "ANY", e2 = cl),  
               function(e1,e2) as(e1,"Matrix") & e2)  
     for(c2 in c("denseMatrix", "Matrix")) {  
         setMethod("&", signature(e1 = cl, e2 = c2),  
                   function(e1,e2) e1 & Diagonal(x = diag(e2)))  
         setMethod("&", signature(e1 = c2, e2 = cl),  
                   function(e1,e2) Diagonal(x = diag(e1)) & e2)  
955      }      }
956  }  }
957    
958    
959  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"
960  ### ----------  any, all: separately here  ### ----------   the last 4: separately here
961  for(cl in diCls) {  for(cl in diCls) {
962  setMethod("any", cl,  setMethod("any", cl,
963            function (x, ..., na.rm) {            function (x, ..., na.rm) {
964                if(any(x@Dim == 0)) FALSE                if(any(x@Dim == 0)) FALSE
965                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)
966            })            })
967  setMethod("all",  cl, function (x, ..., na.rm) any(x@Dim == 0))  setMethod("all",  cl, function (x, ..., na.rm) {
968  setMethod("prod", cl, function (x, ..., na.rm) as.numeric(any(x@Dim == 0)))      n <- x@Dim[1]
969        if(n >= 2) FALSE
970        else if(n == 0 || x@diag == "U") TRUE
971        else all(x@x, ..., na.rm = na.rm)
972    })
973    setMethod("prod", cl, function (x, ..., na.rm) {
974        n <- x@Dim[1]
975        if(n >= 2) 0
976        else if(n == 0 || x@diag == "U") 1
977        else ## n == 1, diag = "N" :
978            prod(x@x, ..., na.rm = na.rm)
979    })
980    
981  setMethod("sum", cl,  setMethod("sum", cl,
982            function(x, ..., na.rm) {            function(x, ..., na.rm) {

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

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