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 2912, Sat Sep 14 17:09:49 2013 UTC revision 3072, Fri Mar 27 15:10:48 2015 UTC
# Line 198  Line 198 
198      as(.bdiag(alis), "CsparseMatrix")      as(.bdiag(alis), "CsparseMatrix")
199  }  }
200    
201    setMethod("tril", "diagonalMatrix", function(x, k = 0, ...)
202        if(k >= 0) x else .setZero(x, paste0(.M.kind(x), "tCMatrix")))
203    
204    setMethod("triu", "diagonalMatrix", function(x, k = 0, ...)
205        if(k <= 0) x else  .setZero(x, paste0(.M.kind(x), "tCMatrix")))
206    
207    
208    
209  .diag2tT <- function(from, uplo = "U", kind = .M.kind(from)) {  .diag2tT <- function(from, uplo = "U", kind = .M.kind(from)) {
210      ## to triangular Tsparse      ## to triangular Tsparse
# Line 227  Line 234 
234                  n))                  n))
235  }  }
236    
237  ## diagonal -> triangular,  upper / lower depending on "partner":  ## diagonal -> triangular,  upper / lower depending on "partner" 'x':
238  diag2tT.u <- function(d, x, kind = .M.kind(d))  diag2tT.u <- function(d, x, kind = .M.kind(d))
239      .diag2tT(d, uplo = if(is(x,"triangularMatrix")) x@uplo else "U", kind)      .diag2tT(d, uplo = if(is(x,"triangularMatrix")) x@uplo else "U", kind)
240    
# Line 251  Line 258 
258  ## "hack"  instead of signature  x = "diagonalMatrix" :  ## "hack"  instead of signature  x = "diagonalMatrix" :
259  ##  ##
260  ## ddi*:  ## ddi*:
261  diag2tT <- function(from) .diag2tT(from, "U", "d")  di2tT <- function(from) .diag2tT(from, "U", "d")
262  setAs("ddiMatrix", "triangularMatrix", diag2tT)  setAs("ddiMatrix", "triangularMatrix", di2tT)
263  ##_no_longer_ setAs("ddiMatrix", "sparseMatrix", diag2tT)  ##_no_longer_ setAs("ddiMatrix", "sparseMatrix", di2tT)
264  ## needed too (otherwise <dense> -> Tsparse is taken):  ## needed too (otherwise <dense> -> Tsparse is taken):
265  setAs("ddiMatrix", "TsparseMatrix", diag2tT)  setAs("ddiMatrix", "TsparseMatrix", di2tT)
266  setAs("ddiMatrix", "dsparseMatrix", diag2tT)  setAs("ddiMatrix", "dsparseMatrix", di2tT)
267  setAs("ddiMatrix", "CsparseMatrix",  setAs("ddiMatrix", "CsparseMatrix",
268        function(from) as(.diag2tT(from, "U", "d"), "CsparseMatrix"))        function(from) as(.diag2tT(from, "U", "d"), "CsparseMatrix"))
269  setAs("ddiMatrix", "symmetricMatrix",  setAs("ddiMatrix", "symmetricMatrix", function(from) .diag2sT(from, "U", "d"))
       function(from) .diag2sT(from, "U", "d"))  
270  ##  ##
271  ## ldi*:  ## ldi*:
272  diag2tT <- function(from) .diag2tT(from, "U", "l")  di2tT <- function(from) .diag2tT(from, "U", "l")
273  setAs("ldiMatrix", "triangularMatrix", diag2tT)  setAs("ldiMatrix", "triangularMatrix", di2tT)
274  ##_no_longer_ setAs("ldiMatrix", "sparseMatrix", diag2tT)  ##_no_longer_ setAs("ldiMatrix", "sparseMatrix", di2tT)
275  ## needed too (otherwise <dense> -> Tsparse is taken):  ## needed too (otherwise <dense> -> Tsparse is taken):
276  setAs("ldiMatrix", "TsparseMatrix", diag2tT)  setAs("ldiMatrix", "TsparseMatrix", di2tT)
277  setAs("ldiMatrix", "lsparseMatrix", diag2tT)  setAs("ldiMatrix", "lsparseMatrix", di2tT)
278  setAs("ldiMatrix", "CsparseMatrix",  setAs("ldiMatrix", "CsparseMatrix",
279        function(from) as(.diag2tT(from, "U", "l"), "CsparseMatrix"))        function(from) as(.diag2tT(from, "U", "l"), "CsparseMatrix"))
280  setAs("ldiMatrix", "symmetricMatrix",  setAs("ldiMatrix", "symmetricMatrix", function(from) .diag2sT(from, "U", "l"))
281        function(from) .diag2sT(from, "U", "l"))  rm(di2tT)
   
282    
283  setAs("diagonalMatrix", "nMatrix",  setAs("diagonalMatrix", "nMatrix",
284        function(from) {        function(from) {
           n <- from@Dim[1]  
285            i <- if(from@diag == "U") integer(0) else which(isN0(from@x)) - 1L            i <- if(from@diag == "U") integer(0) else which(isN0(from@x)) - 1L
286            new("ntTMatrix", i = i, j = i, diag = from@diag,            new("ntTMatrix", i = i, j = i, diag = from@diag,
287                Dim = from@Dim, Dimnames = from@Dimnames)                Dim = from@Dim, Dimnames = from@Dimnames)
# Line 285  Line 289 
289    
290  setAs("diagonalMatrix", "nsparseMatrix", function(from) as(from, "nMatrix"))  setAs("diagonalMatrix", "nsparseMatrix", function(from) as(from, "nMatrix"))
291    
292  ## Cheap fast substitute for diag() which *does* preserve the mode of x :  ##' A version of diag(x,n) which *does* preserve the mode of x, where diag() "fails"
293  mkDiag <- function(x, n) {  mkDiag <- function(x, n) {
294      y <- matrix(as0(mod=mode(x)), n,n)      y <- matrix(as0(mod=mode(x)), n,n)
295      if (n > 0) y[1L + 0:(n - 1L) * (n + 1)] <- x      if (n > 0) y[1L + 0:(n - 1L) * (n + 1)] <- x
296      y      y
297  }  }
298    ## NB: diag(x,n) is really faster for n >= 20, and even more for large n
299    ## --> using diag() where possible, ==> .ddi2mat()
300    
301  setAs("diagonalMatrix", "matrix",  .diag2mat <- function(from)
302        function(from) {      ## want "ldiMatrix" -> <logical> "matrix"  (but integer -> <double> for now)
303            ## want "ldiMatrix" -> <logical> "matrix" :      mkDiag(if(from@diag == "U") as1(from@x) else from@x, n = from@Dim[1])
304            mkDiag(if(from@diag == "U") as1(from@x) else from@x,  
305                   n = from@Dim[1])  .ddi2mat <- function(from)
306        })      base::diag(if(from@diag == "U") as1(from@x) else from@x, nrow = from@Dim[1])
307    
308    setAs("ddiMatrix", "matrix", .ddi2mat)
309    ## the non-ddi diagonalMatrix -- only "ldiMatrix" currently:
310    setAs("diagonalMatrix", "matrix", .diag2mat)
311    
312  setMethod("as.vector", signature(x = "diagonalMatrix", mode="missing"),  setMethod("as.vector", signature(x = "diagonalMatrix", mode="missing"),
313            function(x, mode) {            function(x, mode) {
# Line 516  Line 526 
526            function(x) { x@Dimnames <- x@Dimnames[2:1] ; x })            function(x) { x@Dimnames <- x@Dimnames[2:1] ; x })
527    
528  setMethod("isDiagonal",   "diagonalMatrix", function(object) TRUE)  setMethod("isDiagonal",   "diagonalMatrix", function(object) TRUE)
529  setMethod("isTriangular", "diagonalMatrix", function(object, ...) TRUE)  setMethod("isTriangular", "diagonalMatrix", function(object, upper=NA, ...) TRUE)
530  setMethod("isSymmetric",  "diagonalMatrix", function(object, ...) TRUE)  setMethod("isSymmetric",  "diagonalMatrix", function(object, ...) TRUE)
531    
532  setMethod("symmpart", signature(x = "diagonalMatrix"), function(x) x)  setMethod("symmpart", signature(x = "diagonalMatrix"), function(x) x)
533  setMethod("skewpart", signature(x = "diagonalMatrix"), setZero)  setMethod("skewpart", signature(x = "diagonalMatrix"), function(x) .setZero(x))
534    
535  setMethod("chol", signature(x = "ddiMatrix"),  setMethod("chol", signature(x = "ddiMatrix"),
536            function(x, pivot, ...) {            function(x, pivot, ...) {
# Line 567  Line 577 
577          return(y)          return(y)
578  }  }
579    
580    ##' Boolean Algebra/Arithmetic Product of Diagonal Matrices
581    ##'  %&%
582    diagdiagprodBool <- function(x, y) {
583        dimCheck(x,y)
584        if(x@diag != "U") {
585            if(!is.logical(x@x)) x <- as(x, "lMatrix")
586            if(y@diag != "U") {
587                nx <- x@x & y@x
588                x@x <- as.logical(nx)
589            }
590            x
591        } else { ## x is unit diagonal: return y
592            if(!is.logical(y@x)) y <- as(y, "lMatrix")
593            y
594        }
595    }
596    
597  setMethod("%*%", signature(x = "diagonalMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "diagonalMatrix"),
598            diagdiagprod, valueClass = "ddiMatrix")            diagdiagprod, valueClass = "ddiMatrix")
599    
600  formals(diagdiagprod) <- alist(x=, y=x)  setMethod("%&%", signature(x = "diagonalMatrix", y = "diagonalMatrix"),
601              diagdiagprodBool, valueClass = "ldiMatrix")# do *not* have "ndiMatrix" !
602    
603    ##' Both Numeric or Boolean Algebra/Arithmetic Product of Diagonal Matrices
604    diagdiagprodFlexi <- function(x, y=NULL, boolArith = NA, ...)
605    {
606        dimCheck(x,y)
607        bool <- isTRUE(boolArith)
608        if(x@diag != "U") {
609            if(bool && !is.logical(x@x)) x <- as(x, "lMatrix")
610            if(y@diag != "U") {
611                if(bool) {
612                    nx <- x@x & y@x
613                    x@x <- as.logical(nx)
614                } else { ## boolArith is NA or FALSE: ==> numeric, as have *no* "diagMatrix" patter[n]:
615                    nx <- x@x * y@x
616                    if(is.numeric(nx) && !is.numeric(x@x))
617                        x <- as(x, "dMatrix")
618                    x@x <- as.numeric(nx)
619                }
620            }
621            x
622        } else { ## x is unit diagonal: return y
623            if(bool && !is.logical(y@x)) y <- as(y, "lMatrix")
624            y
625        }
626    }
627  setMethod("crossprod", signature(x = "diagonalMatrix", y = "diagonalMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "diagonalMatrix"),
628            diagdiagprod, valueClass = "ddiMatrix")            diagdiagprodFlexi)
629  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "diagonalMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "diagonalMatrix"),
630            diagdiagprod, valueClass = "ddiMatrix")            diagdiagprodFlexi)
 setMethod("crossprod", signature(x = "diagonalMatrix", y = "missing"),  
           diagdiagprod, valueClass = "ddiMatrix")  
 setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "missing"),  
           diagdiagprod, valueClass = "ddiMatrix")  
631    
632    ##' crossprod(x) := x'x
633    diagprod <- function(x, y = NULL, boolArith = NA, ...) {
634        bool <- isTRUE(boolArith)
635        if(bool && !is.logical(x@x)) x <- as(x, "lMatrix")
636        if(x@diag != "U") {
637            if(bool) {
638                nx <- x@x & y@x
639                x@x <- as.logical(nx)
640            } else { ## boolArith is NA or FALSE: ==> numeric, as have *no* "diagMatrix" patter[n]:
641                nx <- x@x * x@x
642                if(is.numeric(nx) && !is.numeric(x@x))
643                    x <- as(x, "dMatrix")
644                x@x <- as.numeric(nx)
645            }
646        }
647        x
648    }
649    setMethod( "crossprod", signature(x = "diagonalMatrix", y = "missing"), diagprod)
650    setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "missing"), diagprod)
651    
652    
653    ## analogous to matdiagprod() below:
654  diagmatprod <- function(x, y) {  diagmatprod <- function(x, y) {
655      ## x is diagonalMatrix      ## x is diagonalMatrix
656      dx <- dim(x)      if(x@Dim[2] != nrow(y)) stop("non-matching dimensions")
657      dy <- dim(y)      Matrix(if(x@diag == "U") y else x@x * y)
658      if(dx[2] != dy[1]) stop("non-matching dimensions")  }
659      as(if(x@diag == "U") y else x@x * y, "Matrix")  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"), diagmatprod)
660  }  
661  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),  formals(diagmatprod) <- alist(x=, y=NULL, boolArith = NA, ...=)
662            diagmatprod)  setMethod("crossprod",  signature(x = "diagonalMatrix", y = "matrix"), diagmatprod)
 ## sneaky .. :  
 formals(diagmatprod) <- alist(x=, y=NULL)  
 setMethod("crossprod", signature(x = "diagonalMatrix", y = "matrix"),  
           diagmatprod)  
663    
664  diagGeprod <- function(x, y) {  diagGeprod <- function(x, y) {
665      dx <- dim(x)      if(x@Dim[2] != y@Dim[1]) stop("non-matching dimensions")
     dy <- dim(y)  
     if(dx[2] != dy[1]) stop("non-matching dimensions")  
666      if(x@diag != "U")      if(x@diag != "U")
667          y@x <- x@x * y@x          y@x <- x@x * y@x
668      y      y
669  }  }
670  setMethod("%*%", signature(x= "diagonalMatrix", y= "dgeMatrix"), diagGeprod)  setMethod("%*%", signature(x= "diagonalMatrix", y= "dgeMatrix"), diagGeprod)
671  setMethod("%*%", signature(x= "diagonalMatrix", y= "lgeMatrix"), diagGeprod)  setMethod("%*%", signature(x= "diagonalMatrix", y= "lgeMatrix"), diagGeprod)
 formals(diagGeprod) <- alist(x=, y=NULL)  
 setMethod("crossprod", signature(x = "diagonalMatrix", y = "dgeMatrix"),  
           diagGeprod, valueClass = "dgeMatrix")  
 setMethod("crossprod", signature(x = "diagonalMatrix", y = "lgeMatrix"),  
           diagGeprod)  
672    
673    diagGeprodBool <- function(x, y) {
674        if(x@Dim[2] != y@Dim[1]) stop("non-matching dimensions")
675        if(!is.logical(y@x)) y <- as(y, "lMatrix")
676        if(x@diag != "U")
677            y@x <- x@x & y@x
678        y
679    }
680    setMethod("%&%", signature(x= "diagonalMatrix", y= "geMatrix"), diagGeprodBool)
681    
682    diagGeprod2 <- function(x, y=NULL, boolArith = NA, ...) {
683        if(x@Dim[2] != y@Dim[1]) stop("non-matching dimensions")
684        bool <- isTRUE(boolArith)
685        if(bool && !is.logical(y@x)) y <- as(y, "lMatrix")
686        if(x@diag != "U")
687            y@x <- if(bool) x@x & y@x else x@x * y@x
688        y
689    }
690    setMethod("crossprod", signature(x = "diagonalMatrix", y = "dgeMatrix"), diagGeprod2)
691    setMethod("crossprod", signature(x = "diagonalMatrix", y = "lgeMatrix"), diagGeprod2)
692    
693    
694    ## analogous to diagmatprod() above:
695  matdiagprod <- function(x, y) {  matdiagprod <- function(x, y) {
696      dx <- dim(x)      dx <- dim(x)
697      dy <- dim(y)      if(dx[2] != y@Dim[1]) stop("non-matching dimensions")
     if(dx[2] != dy[1]) stop("non-matching dimensions")  
698      Matrix(if(y@diag == "U") x else x * rep(y@x, each = dx[1]))      Matrix(if(y@diag == "U") x else x * rep(y@x, each = dx[1]))
699  }  }
700  setMethod("%*%", signature(x = "matrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "matrix", y = "diagonalMatrix"), matdiagprod)
           matdiagprod)  
 formals(matdiagprod) <- alist(x=, y=NULL)  
 setMethod("tcrossprod", signature(x = "matrix", y = "diagonalMatrix"),  
           matdiagprod)  
701    
702  gediagprod <- function(x, y) {  gediagprod <- function(x, y) {
703      dx <- dim(x)      dx <- dim(x)
704      dy <- dim(y)      if(dx[2] != y@Dim[1]) stop("non-matching dimensions")
     if(dx[2] != dy[1]) stop("non-matching dimensions")  
705      if(y@diag == "N")      if(y@diag == "N")
706          x@x <- x@x * rep(y@x, each = dx[1])          x@x <- x@x * rep(y@x, each = dx[1])
707      x      x
708  }  }
709  setMethod("%*%", signature(x= "dgeMatrix", y= "diagonalMatrix"), gediagprod)  setMethod("%*%", signature(x= "dgeMatrix", y= "diagonalMatrix"), gediagprod)
710  setMethod("%*%", signature(x= "lgeMatrix", y= "diagonalMatrix"), gediagprod)  setMethod("%*%", signature(x= "lgeMatrix", y= "diagonalMatrix"), gediagprod)
711  formals(gediagprod) <- alist(x=, y=NULL)  
712  setMethod("tcrossprod", signature(x = "dgeMatrix", y = "diagonalMatrix"),  gediagprodBool <- function(x, y) {
713            gediagprod)      dx <- dim(x)
714  setMethod("tcrossprod", signature(x = "lgeMatrix", y = "diagonalMatrix"),      if(dx[2] != y@Dim[1]) stop("non-matching dimensions")
715            gediagprod)      if(!is.logical(x@x)) x <- as(x, "lMatrix")
716        if(y@diag == "N")
717            x@x <- x@x & rep(y@x, each = dx[1])
718        x
719    }
720    setMethod("%&%", signature(x= "geMatrix", y= "diagonalMatrix"), gediagprodBool)
721    
722    setMethod("tcrossprod",signature(x = "matrix", y = "diagonalMatrix"),
723              function(x, y=NULL, boolArith = NA, ...) {
724                  dx <- dim(x)
725                  if(dx[2] != y@Dim[1]) stop("non-matching dimensions")
726                  bool <- isTRUE(boolArith)
727                  if(bool && !is.logical(y@x)) y <- as(y, "lMatrix")
728                  Matrix(if(y@diag == "U") x else
729                         if(bool) x & rep(y@x, each = dx[1])
730                         else     x * rep(y@x, each = dx[1]))
731              })
732    
733    setMethod("crossprod", signature(x = "matrix", y = "diagonalMatrix"),
734              function(x, y=NULL, boolArith = NA, ...) {
735                  dx <- dim(x)
736                  if(dx[1] != y@Dim[1]) stop("non-matching dimensions")
737                  bool <- isTRUE(boolArith)
738                  if(bool && !is.logical(y@x)) y <- as(y, "lMatrix")
739                  Matrix(if(bool) t(rep.int(y@x, dx[2]) & x)
740                             else t(rep.int(y@x, dx[2]) * x))
741              })
742    
743    
744    gediagprod2 <- function(x, y=NULL, boolArith = NA, ...) {
745        dx <- dim(x)
746        if(dx[2] != y@Dim[1]) stop("non-matching dimensions")
747        bool <- isTRUE(boolArith)
748        if(bool && !is.logical(x@x)) x <- as(x, "lMatrix")
749        if(y@diag == "N")
750            x@x <- if(bool) x@x & rep(y@x, each = dx[1])
751                   else     x@x * rep(y@x, each = dx[1])
752        x
753    }
754    setMethod("tcrossprod", signature(x = "dgeMatrix", y = "diagonalMatrix"), gediagprod2)
755    setMethod("tcrossprod", signature(x = "lgeMatrix", y = "diagonalMatrix"), gediagprod2)
756    
757    
758  ## crossprod {more of these}  ## crossprod {more of these}
759    
# Line 657  Line 773 
773  ##' @param x CsparseMatrix  ##' @param x CsparseMatrix
774  ##' @param y diagonalMatrix  ##' @param y diagonalMatrix
775  ##' @return x %*% y  ##' @return x %*% y
776  Cspdiagprod <- function(x, y) {  Cspdiagprod <- function(x, y, boolArith = NA, ...) {
777      dx <- dim(x <- .Call(Csparse_diagU2N, x))      if((m <- ncol(x)) != y@Dim[1]) stop("non-matching dimensions")
     dy <- dim(y)  
     if(dx[2] != dy[1]) stop("non-matching dimensions")  
778      if(y@diag == "N") { ## otherwise: y == Diagonal(n) : multiplication is identity      if(y@diag == "N") { ## otherwise: y == Diagonal(n) : multiplication is identity
779          if(!all(y@x[1L] == y@x[-1L]) && is(x, "symmetricMatrix"))          x <- .Call(Csparse_diagU2N, x)
780            cx <- getClass(class(x))
781            if(!all(y@x[1L] == y@x[-1L]) && extends(cx, "symmetricMatrix"))
782              x <- as(x, "generalMatrix")              x <- as(x, "generalMatrix")
783          ind <- rep.int(seq_len(dx[2]), x@p[-1] - x@p[-dx[2]-1L])          ind <- rep.int(seq_len(m), x@p[-1] - x@p[-m-1L])
784            if(isTRUE(boolArith)) {
785                if(extends(cx, "nMatrix")) x <- as(x, "lMatrix") # so, has y@x
786                x@x <- x@x & y@x[x@i + 1L]
787                if(!extends(cx, "diagonalMatrix")) x <- as(drop0(x), "nMatrix")
788            } else {
789                if(!extends(cx, "dMatrix")) x <- as(x, "dMatrix") # <- FIXME if we have zMatrix
790          x@x <- x@x * y@x[ind]          x@x <- x@x * y@x[ind]
791          if(is(x, "compMatrix") && length(xf <- x@factors)) {          }
792            if(.hasSlot(x, "factors") && length(x@factors)) {# drop cashed ones
793              ## instead of dropping all factors, be smart about some              ## instead of dropping all factors, be smart about some
794              ## TODO ......              ## TODO ......
795              x@factors <- list()              x@factors <- list()
# Line 678  Line 801 
801  ##' @param x diagonalMatrix  ##' @param x diagonalMatrix
802  ##' @param y CsparseMatrix  ##' @param y CsparseMatrix
803  ##' @return x %*% y  ##' @return x %*% y
804  diagCspprod <- function(x, y) {  diagCspprod <- function(x, y, boolArith = NA, ...) {
805      dx <- dim(x)      if(x@Dim[2] != y@Dim[1]) stop("non-matching dimensions")
     dy <- dim(y <- .Call(Csparse_diagU2N, y))  
     if(dx[2] != dy[1]) stop("non-matching dimensions")  
806      if(x@diag == "N") {      if(x@diag == "N") {
807          if(!all(x@x[1L] == x@x[-1L]) && is(y, "symmetricMatrix"))          y <- .Call(Csparse_diagU2N, y)
808            cy <- getClass(class(y))
809            if(!all(x@x[1L] == x@x[-1L]) && extends(cy, "symmetricMatrix"))
810              y <- as(y, "generalMatrix")              y <- as(y, "generalMatrix")
811            if(isTRUE(boolArith)) {
812                if(extends(cy, "nMatrix")) y <- as(y, "lMatrix") # so, has y@x
813                y@x <- y@x & x@x[y@i + 1L]
814                if(!extends(cy, "diagonalMatrix")) y <- as(drop0(y), "nMatrix")
815            } else {
816                if(!extends(cy, "dMatrix")) y <- as(y, "dMatrix") # <- FIXME if we have zMatrix
817          y@x <- y@x * x@x[y@i + 1L]          y@x <- y@x * x@x[y@i + 1L]
         if(is(y, "compMatrix") && length(yf <- y@factors)) {  
             ## TODO  
             if(FALSE) { ## instead of dropping all factors, be smart about some  
                 keep <- character()  
                 if(any(iLU <- names(yf) == "LU")) {  
                     keep <- "LU"  
818                  }                  }
819                  y@factors <- yf[keep]          if(.hasSlot(y, "factors") && length(y@factors)) {
820              } else y@factors <- list() ## for now       ## if(.hasSlot(y, "factors") && length(yf <- y@factors)) { ## -- TODO? --
821                ## instead of dropping all factors, be smart about some
822                ## keep <- character()
823                ## if(any(names(yf) == "LU")) { ## <- not easy: y = P'LUQ,  x y = xP'LUQ => LU ???
824                ##     keep <- "LU"
825                ## }
826                ## y@factors <- yf[keep]
827                y@factors <- list()
828          }          }
829      }      }
830      y      y
831  }  }
832    
833    ## + 'boolArith' argument  { ==> .local() is used in any case; keep formals simple :}
834  setMethod("crossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
835            function(x, y = NULL) diagCspprod(x, y))            function(x, y=NULL, boolArith=NA, ...) diagCspprod(x, y, boolArith=boolArith))
836    
837  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
838            function(x, y = NULL) diagCspprod(x, as(y, "CsparseMatrix")))            function(x, y=NULL, boolArith=NA, ...)
839                  diagCspprod(x, as(y, "CsparseMatrix"), boolArith=boolArith))
840    
841  ## Prefer calling diagCspprod to Cspdiagprod if going to transpose anyway  ## Prefer calling diagCspprod to Cspdiagprod if going to transpose anyway
842  ##  x'y == (y'x)'  ##  x'y == (y'x)'
843  setMethod("crossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),  setMethod("crossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
844            function(x, y = NULL) t(diagCspprod(y, x)))            function(x, y=NULL, boolArith=NA, ...) t(diagCspprod(y, x, boolArith=boolArith)))
845    
846  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
847            function(x, y = NULL) t(diagCspprod(y, as(x, "Csparsematrix"))))            function(x, y=NULL, boolArith=NA, ...) t(diagCspprod(y, as(x, "Csparsematrix"), boolArith=boolArith)))
848    
849  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
850            function(x, y = NULL) diagCspprod(x, t(y)))            function(x, y=NULL, boolArith=NA, ...) diagCspprod(x, t(y), boolArith=boolArith))
851    
852  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
853            function(x, y = NULL) diagCspprod(x, t(as(y, "CsparseMatrix"))))            function(x, y=NULL, boolArith=NA, ...) diagCspprod(x, t(as(y, "CsparseMatrix")), boolArith=boolArith))
854    
855  setMethod("tcrossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),  setMethod("tcrossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
856            function(x, y = NULL) Cspdiagprod(x, y))            function(x, y=NULL, boolArith=NA, ...) Cspdiagprod(x, y, boolArith=boolArith))
857    
858  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
859            function(x, y = NULL) Cspdiagprod(as(x, "CsparseMatrix"), y))            function(x, y=NULL, boolArith=NA, ...) Cspdiagprod(as(x, "CsparseMatrix"), y, boolArith=boolArith))
860    
861  setMethod("%*%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
862            function(x, y) diagCspprod(x, y))            function(x, y) diagCspprod(x, y, boolArith=NA))
863    setMethod("%&%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
864              function(x, y) diagCspprod(x, y, boolArith=TRUE))
865    
866  ## instead of "sparseMatrix", use: [RT]sparse.. ("closer" in method dispatch)  ## instead of "sparseMatrix", use: [RT]sparse.. ("closer" in method dispatch)
867  for(cl in c("TsparseMatrix", "RsparseMatrix")) {  for(cl in c("TsparseMatrix", "RsparseMatrix")) {
868    
869  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
870            function(x, y) diagCspprod(as(x, "CsparseMatrix"), y))            function(x, y) diagCspprod(as(x, "CsparseMatrix"), y, boolArith=NA))
871    
872  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
873            function(x, y) Cspdiagprod(as(x, "CsparseMatrix"), y))            function(x, y) Cspdiagprod(as(x, "CsparseMatrix"), y, boolArith=NA))
874    
875    setMethod("%&%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
876              function(x, y) diagCspprod(as(x, "CsparseMatrix"), y, boolArith=TRUE))
877    
878    setMethod("%&%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
879              function(x, y) Cspdiagprod(as(x, "CsparseMatrix"), y, boolArith=TRUE))
880    
881  }  }
882    
883  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
884            function(x, y) Cspdiagprod(x, y))            function(x, y) Cspdiagprod(x, y, boolArith=NA))
885    setMethod("%&%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
886              function(x, y) Cspdiagprod(x, y, boolArith=TRUE))
887    
888  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*
889  ##       do indeed work by going through sparse (and *not* ddense)!  ##       do indeed work by going through sparse (and *not* ddense)!
# Line 833  Line 976 
976      ## result must be triangular      ## result must be triangular
977      r <- callGeneric(d1 <- .diag.x(e1), diag(e2)) # error if not "compatible"      r <- callGeneric(d1 <- .diag.x(e1), diag(e2)) # error if not "compatible"
978      ## Check what happens with non-diagonals, i.e. (0 o 0), (FALSE o 0), ...:      ## Check what happens with non-diagonals, i.e. (0 o 0), (FALSE o 0), ...:
979      e1.0 <- if(.n1 <- is.numeric(d1   )) 0 else FALSE      e1.0 <- if(is.numeric(d1)) 0 else FALSE
980      r00 <- callGeneric(e1.0, if(.n2 <- is.numeric(e2[0])) 0 else FALSE)      r00 <- callGeneric(e1.0, if(.n2 <- is.numeric(e2[0])) 0 else FALSE)
981      if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* triangular      if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* triangular
982          diag(e2) <- r          diag(e2) <- r
# Line 918  Line 1061 
1061            function(e1,e2) {            function(e1,e2) {
1062                n <- e1@Dim[1]                n <- e1@Dim[1]
1063                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
1064                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
1065                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
1066                    if(e1@diag == "U") {                    if(e1@diag == "U") {
1067                        if(any((r <- callGeneric(1, e2)) != 1)) {                        if(any((r <- callGeneric(1, e2)) != 1)) {
# Line 940  Line 1083 
1083            function(e1,e2) {            function(e1,e2) {
1084                n <- e2@Dim[1]                n <- e2@Dim[1]
1085                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
1086                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
1087                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
1088                    if(e2@diag == "U") {                    if(e2@diag == "U") {
1089                        if(any((r <- callGeneric(e1, 1)) != 1)) {                        if(any((r <- callGeneric(e1, 1)) != 1)) {
# Line 963  Line 1106 
1106            function(e1,e2) {            function(e1,e2) {
1107                n <- e1@Dim[1]                n <- e1@Dim[1]
1108                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
1109                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
1110                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
1111                    E <- copyClass(e1, "ddiMatrix", c("diag", "Dim", "Dimnames"))#FIXME: if ok, check=FALSE                    E <- copyClass(e1, "ddiMatrix", c("diag", "Dim", "Dimnames"), check=FALSE)
                   ## E <- copyClass(e1, "ddiMatrix", check=FALSE)  
1112                    ## storage.mode(E@x) <- "double"                    ## storage.mode(E@x) <- "double"
1113                    if(e1@diag == "U") {                    if(e1@diag == "U") {
1114                        if(any((r <- callGeneric(1, e2)) != 1)) {                        if(any((r <- callGeneric(1, e2)) != 1)) {
# Line 988  Line 1130 
1130            function(e1,e2) {            function(e1,e2) {
1131                n <- e2@Dim[1]                n <- e2@Dim[1]
1132                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
1133                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
1134                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
1135                    E <- copyClass(e2, "ddiMatrix", c("diag", "Dim", "Dimnames"))#FIXME: if ok, check=FALSE                    E <- copyClass(e2, "ddiMatrix", c("diag", "Dim", "Dimnames"), check=FALSE)
                   ## E <- copyClass(e2, "ddiMatrix", check=FALSE)  
1136                    ## storage.mode(E@x) <- "double"                    ## storage.mode(E@x) <- "double"
1137                    if(e2@diag == "U") {                    if(e2@diag == "U") {
1138                        if(any((r <- callGeneric(e1, 1)) != 1)) {                        if(any((r <- callGeneric(e1, 1)) != 1)) {
# Line 1021  Line 1162 
1162            function(e1,e2) {            function(e1,e2) {
1163                n <- e1@Dim[1]                n <- e1@Dim[1]
1164                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
1165                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
1166                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
1167                    E <- copyClass(e1, "ldiMatrix", c("diag", "Dim", "Dimnames"))#FIXME: if ok, check=FALSE                    E <- copyClass(e1, "ldiMatrix", c("diag", "Dim", "Dimnames"), check=FALSE)
                   ## E <- copyClass(e1, "ldiMatrix", check=FALSE)  
1168                    ## storage.mode(E@x) <- "logical"                    ## storage.mode(E@x) <- "logical"
1169                    if(e1@diag == "U") {                    if(e1@diag == "U") {
1170                        if(any((r <- callGeneric(1, e2)) != 1)) {                        if(any((r <- callGeneric(1, e2)) != 1)) {
# Line 1047  Line 1187 
1187            function(e1,e2) {            function(e1,e2) {
1188                n <- e1@Dim[1]                n <- e1@Dim[1]
1189                f0 <- callGeneric(FALSE, e2)                f0 <- callGeneric(FALSE, e2)
1190                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
1191                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
1192    
1193                    if(e1@diag == "U") {                    if(e1@diag == "U") {
# Line 1081  Line 1221 
1221  }  }
1222    
1223  ## Direct subclasses of "denseMatrix": currently ddenseMatrix, ldense... :  ## Direct subclasses of "denseMatrix": currently ddenseMatrix, ldense... :
1224    if(FALSE) # now also contains "geMatrix"
1225  dense.subCl <- local({ dM.scl <- getClass("denseMatrix")@subclasses  dense.subCl <- local({ dM.scl <- getClass("denseMatrix")@subclasses
1226                         names(dM.scl)[vapply(dM.scl, slot, 0, "distance") == 1] })                         names(dM.scl)[vapply(dM.scl, slot, 0, "distance") == 1] })
1227    dense.subCl <- paste0(c("d","l","n"), "denseMatrix")
1228  for(DI in diCls) {  for(DI in diCls) {
1229      dMeth <- if(extends(DI, "dMatrix"))      dMeth <- if(extends(DI, "dMatrix"))
1230          function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2)          function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2)
# Line 1102  Line 1244 
1244      }      }
1245  }  }
1246    
1247    ## Group methods "Math", "Math2" in                     --> ./Math.R
1248    
1249  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"
1250  ### ----------   the last 4: separately here  ### ----------   the last 4: separately here
# Line 1181  Line 1324 
1324                }                }
1325            })            })
1326    
1327  rm(arg1, arg2, other, DI, cl, c1, c2,  rm(arg1, arg2, other, DI, Fun, cl, c1, c2,
1328     dense.subCl, diCls)# not used elsewhere     dense.subCl, diCls)# not used elsewhere
1329    
1330  setMethod("summary", signature(object = "diagonalMatrix"),  setMethod("summary", signature(object = "diagonalMatrix"),

Legend:
Removed from v.2912  
changed lines
  Added in v.3072

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