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 2922, Thu Sep 26 16:36:49 2013 UTC revision 3079, Tue Mar 31 15:29:43 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 259  Line 266 
266  setAs("ddiMatrix", "dsparseMatrix", di2tT)  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  di2tT <- function(from) .diag2tT(from, "U", "l")  di2tT <- function(from) .diag2tT(from, "U", "l")
# Line 271  Line 277 
277  setAs("ldiMatrix", "lsparseMatrix", di2tT)  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"))
       function(from) .diag2sT(from, "U", "l"))  
281  rm(di2tT)  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, ...=) ## FIXME boolArith
662            diagmatprod)  diagmatprod2 <- function(x, y=NULL, boolArith = NA, ...) {
663  ## sneaky .. :      ## x is diagonalMatrix
664  formals(diagmatprod) <- alist(x=, y=NULL)      if(x@Dim[2] != nrow(y)) stop("non-matching dimensions")
665  setMethod("crossprod", signature(x = "diagonalMatrix", y = "matrix"),      Matrix(if(x@diag == "U") y else x@x * y)
666            diagmatprod)  }
667    setMethod("crossprod",  signature(x = "diagonalMatrix", y = "matrix"), diagmatprod2)
668    
669  diagGeprod <- function(x, y) {  diagGeprod <- function(x, y) {
670      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")  
671      if(x@diag != "U")      if(x@diag != "U")
672          y@x <- x@x * y@x          y@x <- x@x * y@x
673      y      y
674  }  }
675  setMethod("%*%", signature(x= "diagonalMatrix", y= "dgeMatrix"), diagGeprod)  setMethod("%*%", signature(x= "diagonalMatrix", y= "dgeMatrix"), diagGeprod)
676  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)  
677    
678    diagGeprodBool <- function(x, y) {
679        if(x@Dim[2] != y@Dim[1]) stop("non-matching dimensions")
680        if(!is.logical(y@x)) y <- as(y, "lMatrix")
681        if(x@diag != "U")
682            y@x <- x@x & y@x
683        y
684    }
685    setMethod("%&%", signature(x= "diagonalMatrix", y= "geMatrix"), diagGeprodBool)
686    
687    diagGeprod2 <- function(x, y=NULL, boolArith = NA, ...) {
688        if(x@Dim[2] != y@Dim[1]) stop("non-matching dimensions")
689        bool <- isTRUE(boolArith)
690        if(bool && !is.logical(y@x)) y <- as(y, "lMatrix")
691        if(x@diag != "U")
692            y@x <- if(bool) x@x & y@x else x@x * y@x
693        y
694    }
695    setMethod("crossprod", signature(x = "diagonalMatrix", y = "dgeMatrix"), diagGeprod2)
696    setMethod("crossprod", signature(x = "diagonalMatrix", y = "lgeMatrix"), diagGeprod2)
697    
698    
699    ## analogous to diagmatprod() above:
700  matdiagprod <- function(x, y) {  matdiagprod <- function(x, y) {
701      dx <- dim(x)      dx <- dim(x)
702      dy <- dim(y)      if(dx[2] != y@Dim[1]) stop("non-matching dimensions")
     if(dx[2] != dy[1]) stop("non-matching dimensions")  
703      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]))
704  }  }
705  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)  
706    
707  gediagprod <- function(x, y) {  gediagprod <- function(x, y) {
708      dx <- dim(x)      dx <- dim(x)
709      dy <- dim(y)      if(dx[2] != y@Dim[1]) stop("non-matching dimensions")
     if(dx[2] != dy[1]) stop("non-matching dimensions")  
710      if(y@diag == "N")      if(y@diag == "N")
711          x@x <- x@x * rep(y@x, each = dx[1])          x@x <- x@x * rep(y@x, each = dx[1])
712      x      x
713  }  }
714  setMethod("%*%", signature(x= "dgeMatrix", y= "diagonalMatrix"), gediagprod)  setMethod("%*%", signature(x= "dgeMatrix", y= "diagonalMatrix"), gediagprod)
715  setMethod("%*%", signature(x= "lgeMatrix", y= "diagonalMatrix"), gediagprod)  setMethod("%*%", signature(x= "lgeMatrix", y= "diagonalMatrix"), gediagprod)
716  formals(gediagprod) <- alist(x=, y=NULL)  
717  setMethod("tcrossprod", signature(x = "dgeMatrix", y = "diagonalMatrix"),  gediagprodBool <- function(x, y) {
718            gediagprod)      dx <- dim(x)
719  setMethod("tcrossprod", signature(x = "lgeMatrix", y = "diagonalMatrix"),      if(dx[2] != y@Dim[1]) stop("non-matching dimensions")
720            gediagprod)      if(!is.logical(x@x)) x <- as(x, "lMatrix")
721        if(y@diag == "N")
722            x@x <- x@x & rep(y@x, each = dx[1])
723        x
724    }
725    setMethod("%&%", signature(x= "geMatrix", y= "diagonalMatrix"), gediagprodBool)
726    
727    setMethod("tcrossprod",signature(x = "matrix", y = "diagonalMatrix"),
728              function(x, y=NULL, boolArith = NA, ...) {
729                  dx <- dim(x)
730                  if(dx[2] != y@Dim[1]) stop("non-matching dimensions")
731                  bool <- isTRUE(boolArith)
732                  if(bool && !is.logical(y@x)) y <- as(y, "lMatrix")
733                  Matrix(if(y@diag == "U") x else
734                         if(bool) x & rep(y@x, each = dx[1])
735                         else     x * rep(y@x, each = dx[1]))
736              })
737    
738    setMethod("crossprod", signature(x = "matrix", y = "diagonalMatrix"),
739              function(x, y=NULL, boolArith = NA, ...) {
740                  dx <- dim(x)
741                  if(dx[1] != y@Dim[1]) stop("non-matching dimensions")
742                  bool <- isTRUE(boolArith)
743                  if(bool && !is.logical(y@x)) y <- as(y, "lMatrix")
744                  Matrix(if(bool) t(rep.int(y@x, dx[2]) & x)
745                             else t(rep.int(y@x, dx[2]) * x))
746              })
747    
748    
749    gediagprod2 <- function(x, y=NULL, boolArith = NA, ...) {
750        dx <- dim(x)
751        if(dx[2] != y@Dim[1]) stop("non-matching dimensions")
752        bool <- isTRUE(boolArith)
753        if(bool && !is.logical(x@x)) x <- as(x, "lMatrix")
754        if(y@diag == "N")
755            x@x <- if(bool) x@x & rep(y@x, each = dx[1])
756                   else     x@x * rep(y@x, each = dx[1])
757        x
758    }
759    setMethod("tcrossprod", signature(x = "dgeMatrix", y = "diagonalMatrix"), gediagprod2)
760    setMethod("tcrossprod", signature(x = "lgeMatrix", y = "diagonalMatrix"), gediagprod2)
761    
762    
763  ## crossprod {more of these}  ## crossprod {more of these}
764    
# Line 657  Line 778 
778  ##' @param x CsparseMatrix  ##' @param x CsparseMatrix
779  ##' @param y diagonalMatrix  ##' @param y diagonalMatrix
780  ##' @return x %*% y  ##' @return x %*% y
781  Cspdiagprod <- function(x, y) {  Cspdiagprod <- function(x, y, boolArith = NA, ...) {
782      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")  
783      if(y@diag == "N") { ## otherwise: y == Diagonal(n) : multiplication is identity      if(y@diag == "N") { ## otherwise: y == Diagonal(n) : multiplication is identity
784          if(!all(y@x[1L] == y@x[-1L]) && is(x, "symmetricMatrix"))          x <- .Call(Csparse_diagU2N, x)
785            cx <- getClass(class(x))
786            if(!all(y@x[1L] == y@x[-1L]) && extends(cx, "symmetricMatrix"))
787              x <- as(x, "generalMatrix")              x <- as(x, "generalMatrix")
788          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])
789            if(isTRUE(boolArith)) {
790                if(extends(cx, "nMatrix")) x <- as(x, "lMatrix") # so, has y@x
791                x@x <- r <- x@x & y@x[x@i + 1L]
792                if(!anyNA(r) && !extends(cx, "diagonalMatrix")) x <- as(drop0(x), "nMatrix")
793            } else {
794                if(!extends(cx, "dMatrix")) x <- as(x, "dMatrix") # <- FIXME if we have zMatrix
795          x@x <- x@x * y@x[ind]          x@x <- x@x * y@x[ind]
796          if(is(x, "compMatrix") && length(xf <- x@factors)) {          }
797            if(.hasSlot(x, "factors") && length(x@factors)) {# drop cashed ones
798              ## instead of dropping all factors, be smart about some              ## instead of dropping all factors, be smart about some
799              ## TODO ......              ## TODO ......
800              x@factors <- list()              x@factors <- list()
801          }          }
     }  
802      x      x
803        } else { #  y is unit-diagonal ==> "return x"
804            cx <- getClass(class(x))
805            if(isTRUE(boolArith)) {
806                is.l <- if(extends(cx, "dMatrix")) { ## <- FIXME: extend once we have iMatrix, zMatrix
807                    x <- as(x, "lMatrix"); TRUE } else extends(cx, "lMatrix")
808                if(is.l && !anyNA(x@x)) as(drop0(x), "nMatrix")
809                else if(is.l) x else # defensive:
810                as(x, "lMatrix")
811            } else {
812                ## else boolArith is  NA or FALSE {which are equivalent here, das diagonal = "numLike"}
813                if(extends(cx, "nMatrix") || extends(cx, "lMatrix"))
814                    as(x, "dMatrix") else x
815            }
816        }
817  }  }
818    
819  ##' @param x diagonalMatrix  ##' @param x diagonalMatrix
820  ##' @param y CsparseMatrix  ##' @param y CsparseMatrix
821  ##' @return x %*% y  ##' @return x %*% y
822  diagCspprod <- function(x, y) {  diagCspprod <- function(x, y, boolArith = NA, ...) {
823      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")  
824      if(x@diag == "N") {      if(x@diag == "N") {
825          if(!all(x@x[1L] == x@x[-1L]) && is(y, "symmetricMatrix"))          y <- .Call(Csparse_diagU2N, y)
826            cy <- getClass(class(y))
827            if(!all(x@x[1L] == x@x[-1L]) && extends(cy, "symmetricMatrix"))
828              y <- as(y, "generalMatrix")              y <- as(y, "generalMatrix")
829            if(isTRUE(boolArith)) {
830                if(extends(cy, "nMatrix")) y <- as(y, "lMatrix") # so, has y@x
831                y@x <- r <- y@x & x@x[y@i + 1L]
832                if(!anyNA(r) && !extends(cy, "diagonalMatrix")) y <- as(drop0(y), "nMatrix")
833            } else {
834                if(!extends(cy, "dMatrix")) y <- as(y, "dMatrix") # <- FIXME if we have zMatrix
835          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"  
                 }  
                 y@factors <- yf[keep]  
             } else y@factors <- list() ## for now  
836          }          }
837            if(.hasSlot(y, "factors") && length(y@factors)) {
838         ## if(.hasSlot(y, "factors") && length(yf <- y@factors)) { ## -- TODO? --
839                ## instead of dropping all factors, be smart about some
840                ## keep <- character()
841                ## if(any(names(yf) == "LU")) { ## <- not easy: y = P'LUQ,  x y = xP'LUQ => LU ???
842                ##     keep <- "LU"
843                ## }
844                ## y@factors <- yf[keep]
845                y@factors <- list()
846      }      }
847      y      y
848        } else { ## x @ diag  == "U"
849            cy <- getClass(class(y))
850            if(isTRUE(boolArith)) {
851                is.l <- if(extends(cy, "dMatrix")) { ## <- FIXME: extend once we have iMatrix, zMatrix
852                    y <- as(y, "lMatrix"); TRUE } else extends(cy, "lMatrix")
853                if(is.l && !anyNA(y@x)) as(drop0(y), "nMatrix")
854                else if(is.l) y else # defensive:
855                as(y, "lMatrix")
856            } else {
857                ## else boolArith is  NA or FALSE {which are equivalent here, das diagonal = "numLike"}
858                if(extends(cy, "nMatrix") || extends(cy, "lMatrix"))
859                    as(y, "dMatrix") else y
860            }
861        }
862  }  }
863    
864    ## + 'boolArith' argument  { ==> .local() is used in any case; keep formals simple :}
865  setMethod("crossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
866            function(x, y = NULL) diagCspprod(x, y))            function(x, y=NULL, boolArith=NA, ...) diagCspprod(x, y, boolArith=boolArith))
867    
868  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
869            function(x, y = NULL) diagCspprod(x, as(y, "CsparseMatrix")))            function(x, y=NULL, boolArith=NA, ...)
870                  diagCspprod(x, as(y, "CsparseMatrix"), boolArith=boolArith))
871    
872  ## Prefer calling diagCspprod to Cspdiagprod if going to transpose anyway  ## Prefer calling diagCspprod to Cspdiagprod if going to transpose anyway
873  ##  x'y == (y'x)'  ##  x'y == (y'x)'
874  setMethod("crossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),  setMethod("crossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
875            function(x, y = NULL) t(diagCspprod(y, x)))            function(x, y=NULL, boolArith=NA, ...) t(diagCspprod(y, x, boolArith=boolArith)))
876    
877  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
878            function(x, y = NULL) t(diagCspprod(y, as(x, "Csparsematrix"))))            function(x, y=NULL, boolArith=NA, ...) t(diagCspprod(y, as(x, "Csparsematrix"), boolArith=boolArith)))
879    
880  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
881            function(x, y = NULL) diagCspprod(x, t(y)))            function(x, y=NULL, boolArith=NA, ...) diagCspprod(x, t(y), boolArith=boolArith))
882    
883  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
884            function(x, y = NULL) diagCspprod(x, t(as(y, "CsparseMatrix"))))            function(x, y=NULL, boolArith=NA, ...) diagCspprod(x, t(as(y, "CsparseMatrix")), boolArith=boolArith))
885    
886  setMethod("tcrossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),  setMethod("tcrossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
887            function(x, y = NULL) Cspdiagprod(x, y))            function(x, y=NULL, boolArith=NA, ...) Cspdiagprod(x, y, boolArith=boolArith))
888    
889  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
890            function(x, y = NULL) Cspdiagprod(as(x, "CsparseMatrix"), y))            function(x, y=NULL, boolArith=NA, ...) Cspdiagprod(as(x, "CsparseMatrix"), y, boolArith=boolArith))
891    
892  setMethod("%*%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
893            function(x, y) diagCspprod(x, y))            function(x, y) diagCspprod(x, y, boolArith=NA))
894    setMethod("%&%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
895              function(x, y) diagCspprod(x, y, boolArith=TRUE))
896    
897  ## instead of "sparseMatrix", use: [RT]sparse.. ("closer" in method dispatch)  ## instead of "sparseMatrix", use: [RT]sparse.. ("closer" in method dispatch)
898  for(cl in c("TsparseMatrix", "RsparseMatrix")) {  for(cl in c("TsparseMatrix", "RsparseMatrix")) {
899    
900  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
901            function(x, y) diagCspprod(as(x, "CsparseMatrix"), y))            function(x, y) diagCspprod(as(x, "CsparseMatrix"), y, boolArith=NA))
902    
903  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
904            function(x, y) Cspdiagprod(as(x, "CsparseMatrix"), y))            function(x, y) Cspdiagprod(as(x, "CsparseMatrix"), y, boolArith=NA))
905    
906    setMethod("%&%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
907              function(x, y) diagCspprod(as(x, "CsparseMatrix"), y, boolArith=TRUE))
908    
909    setMethod("%&%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
910              function(x, y) Cspdiagprod(as(x, "CsparseMatrix"), y, boolArith=TRUE))
911    
912  }  }
913    
914  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
915            function(x, y) Cspdiagprod(x, y))            function(x, y) Cspdiagprod(x, y, boolArith=NA))
916    setMethod("%&%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
917              function(x, y) Cspdiagprod(x, y, boolArith=TRUE))
918    
919  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*
920  ##       do indeed work by going through sparse (and *not* ddense)!  ##       do indeed work by going through sparse (and *not* ddense)!
# Line 833  Line 1007 
1007      ## result must be triangular      ## result must be triangular
1008      r <- callGeneric(d1 <- .diag.x(e1), diag(e2)) # error if not "compatible"      r <- callGeneric(d1 <- .diag.x(e1), diag(e2)) # error if not "compatible"
1009      ## 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), ...:
1010      e1.0 <- if(.n1 <- is.numeric(d1   )) 0 else FALSE      e1.0 <- if(is.numeric(d1)) 0 else FALSE
1011      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)
1012      if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* triangular      if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* triangular
1013          diag(e2) <- r          diag(e2) <- r
# Line 918  Line 1092 
1092            function(e1,e2) {            function(e1,e2) {
1093                n <- e1@Dim[1]                n <- e1@Dim[1]
1094                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
1095                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
1096                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
1097                    if(e1@diag == "U") {                    if(e1@diag == "U") {
1098                        if(any((r <- callGeneric(1, e2)) != 1)) {                        if(any((r <- callGeneric(1, e2)) != 1)) {
# Line 940  Line 1114 
1114            function(e1,e2) {            function(e1,e2) {
1115                n <- e2@Dim[1]                n <- e2@Dim[1]
1116                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
1117                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
1118                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
1119                    if(e2@diag == "U") {                    if(e2@diag == "U") {
1120                        if(any((r <- callGeneric(e1, 1)) != 1)) {                        if(any((r <- callGeneric(e1, 1)) != 1)) {
# Line 963  Line 1137 
1137            function(e1,e2) {            function(e1,e2) {
1138                n <- e1@Dim[1]                n <- e1@Dim[1]
1139                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
1140                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
1141                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
1142                    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)  
1143                    ## storage.mode(E@x) <- "double"                    ## storage.mode(E@x) <- "double"
1144                    if(e1@diag == "U") {                    if(e1@diag == "U") {
1145                        if(any((r <- callGeneric(1, e2)) != 1)) {                        if(any((r <- callGeneric(1, e2)) != 1)) {
# Line 988  Line 1161 
1161            function(e1,e2) {            function(e1,e2) {
1162                n <- e2@Dim[1]                n <- e2@Dim[1]
1163                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
1164                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
1165                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
1166                    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)  
1167                    ## storage.mode(E@x) <- "double"                    ## storage.mode(E@x) <- "double"
1168                    if(e2@diag == "U") {                    if(e2@diag == "U") {
1169                        if(any((r <- callGeneric(e1, 1)) != 1)) {                        if(any((r <- callGeneric(e1, 1)) != 1)) {
# Line 1021  Line 1193 
1193            function(e1,e2) {            function(e1,e2) {
1194                n <- e1@Dim[1]                n <- e1@Dim[1]
1195                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
1196                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
1197                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
1198                    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)  
1199                    ## storage.mode(E@x) <- "logical"                    ## storage.mode(E@x) <- "logical"
1200                    if(e1@diag == "U") {                    if(e1@diag == "U") {
1201                        if(any((r <- callGeneric(1, e2)) != 1)) {                        if(any((r <- callGeneric(1, e2)) != 1)) {
# Line 1047  Line 1218 
1218            function(e1,e2) {            function(e1,e2) {
1219                n <- e1@Dim[1]                n <- e1@Dim[1]
1220                f0 <- callGeneric(FALSE, e2)                f0 <- callGeneric(FALSE, e2)
1221                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
1222                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
1223    
1224                    if(e1@diag == "U") {                    if(e1@diag == "U") {
# Line 1081  Line 1252 
1252  }  }
1253    
1254  ## Direct subclasses of "denseMatrix": currently ddenseMatrix, ldense... :  ## Direct subclasses of "denseMatrix": currently ddenseMatrix, ldense... :
1255    if(FALSE) # now also contains "geMatrix"
1256  dense.subCl <- local({ dM.scl <- getClass("denseMatrix")@subclasses  dense.subCl <- local({ dM.scl <- getClass("denseMatrix")@subclasses
1257                         names(dM.scl)[vapply(dM.scl, slot, 0, "distance") == 1] })                         names(dM.scl)[vapply(dM.scl, slot, 0, "distance") == 1] })
1258    dense.subCl <- paste0(c("d","l","n"), "denseMatrix")
1259  for(DI in diCls) {  for(DI in diCls) {
1260      dMeth <- if(extends(DI, "dMatrix"))      dMeth <- if(extends(DI, "dMatrix"))
1261          function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2)          function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2)
# Line 1102  Line 1275 
1275      }      }
1276  }  }
1277    
1278    ## Group methods "Math", "Math2" in                     --> ./Math.R
1279    
1280  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"
1281  ### ----------   the last 4: separately here  ### ----------   the last 4: separately here
# Line 1181  Line 1355 
1355                }                }
1356            })            })
1357    
1358  rm(arg1, arg2, other, DI, cl, c1, c2,  rm(arg1, arg2, other, DI, Fun, cl, c1, c2,
1359     dense.subCl, diCls)# not used elsewhere     dense.subCl, diCls)# not used elsewhere
1360    
1361  setMethod("summary", signature(object = "diagonalMatrix"),  setMethod("summary", signature(object = "diagonalMatrix"),

Legend:
Removed from v.2922  
changed lines
  Added in v.3079

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