--- pkg/R/Matrix.R 2006/08/28 15:35:44 1455 +++ pkg/R/Matrix.R 2007/01/29 20:17:33 1747 @@ -6,6 +6,11 @@ setAs("Matrix", "sparseMatrix", function(from) as_Csparse(from)) setAs("Matrix", "denseMatrix", function(from) as_dense(from)) +## Most of these work; this is a last resort: +setAs(from = "Matrix", to = "matrix", # do *not* call base::as.matrix() here: + function(from) .bail.out.2("coerce", class(from), class(to))) +setAs(from = "matrix", to = "Matrix", function(from) Matrix(from)) + ## ## probably not needed eventually: ## setAs(from = "ddenseMatrix", to = "matrix", ## function(from) { @@ -19,13 +24,8 @@ setMethod("as.array", signature(x = "Matrix"), function(x) as(x, "matrix")) ## head and tail apply to all Matrix objects for which subscripting is allowed: -## if(paste(R.version\$major, R.version\$minor, sep=".") < "2.4") { - setMethod("head", signature(x = "Matrix"), utils:::head.matrix) - setMethod("tail", signature(x = "Matrix"), utils:::tail.matrix) -## } else { # R 2.4.0 and newer -## setMethod("head", signature(x = "Matrix"), utils::head.matrix) -## setMethod("tail", signature(x = "Matrix"), utils::tail.matrix) -## } +setMethod("head", signature(x = "Matrix"), utils::head.matrix) +setMethod("tail", signature(x = "Matrix"), utils::tail.matrix) ## slow "fall back" method {subclasses should have faster ones}: setMethod("as.vector", signature(x = "Matrix", mode = "missing"), @@ -38,34 +38,13 @@ function(x, ...) as.logical(as.vector(x))) -## Note that isSymmetric is *not* exported -## but that "base" has an isSymmetric() S3-generic since R 2.3.0 +## "base" has an isSymmetric() S3-generic since R 2.3.0 setMethod("isSymmetric", signature(object = "symmetricMatrix"), function(object,tol) TRUE) setMethod("isSymmetric", signature(object = "triangularMatrix"), ## TRUE iff diagonal: function(object,tol) isDiagonal(object)) -if(paste(R.version\$major, R.version\$minor, sep=".") < "2.3") - ## need a "matrix" method as in R 2.3 and later - setMethod("isSymmetric", signature(object = "matrix"), - function(object, tol = 100*.Machine\$double.eps, ...) - { - ## pretest: is it square? - d <- dim(object) - if(d[1] != d[2]) return(FALSE) - ## for `broken' all.equal in R <= 2.2.x: - dn <- dimnames(object) - if(!identical(dn[1], dn[2])) return(FALSE) - test <- - if(is.complex(object)) - all.equal.numeric(object, Conj(t(object)), tol = tol, ...) - else # numeric, character, .. - all.equal(object, t(object), tol = tol, ...) - isTRUE(test) - }) - - setMethod("isTriangular", signature(object = "triangularMatrix"), function(object, ...) TRUE) @@ -77,7 +56,12 @@ setMethod("dim", signature(x = "Matrix"), function(x) x@Dim, valueClass = "integer") + +setMethod("length", "Matrix", function(x) prod(dim(x))) + setMethod("dimnames", signature(x = "Matrix"), function(x) x@Dimnames) + + ## not exported but used more than once for "dimnames<-" method : ## -- or do only once for all "Matrix" classes ?? dimnamesGets <- function (x, value) { @@ -96,21 +80,30 @@ setMethod("unname", signature("Matrix", force="missing"), function(obj) { obj@Dimnames <- list(NULL,NULL); obj}) +setMethod("all", signature(x = "Matrix"), + function(x, ..., na.rm) { x <- as(x, "lMatrix"); callGeneric()}) +setMethod("any", signature(x = "Matrix"), + function(x, ..., na.rm) { x <- as(x, "lMatrix"); callGeneric()}) + +setMethod("!", "Matrix", function(e1) !as(e1, "lMatrix")) + + + Matrix <- function (data = NA, nrow = 1, ncol = 1, byrow = FALSE, dimnames = NULL, sparse = NULL, forceCheck = FALSE) { - sparseDefault <- function(m) - prod(dim(m)) > 2*sum(as(m, "matrix") != 0) + sparseDefault <- function(m) prod(dim(m)) > 2*sum(isN0(as(m, "matrix"))) i.M <- is(data, "Matrix") - if(is.null(sparse) && (i.M || is(data, "matrix"))) + + if(is.null(sparse1 <- sparse) && (i.M || is(data, "matrix"))) sparse <- sparseDefault(data) doDN <- TRUE - if (i.M && !forceCheck) { + if (i.M) { sM <- is(data,"sparseMatrix") - if((sparse && sM) || (!sparse && !sM)) + if(!forceCheck && ((sparse && sM) || (!sparse && !sM))) return(data) ## else : convert dense <-> sparse -> at end } @@ -119,8 +112,9 @@ nrow <- ceiling(length(data)/ncol) else if (missing(ncol)) ncol <- ceiling(length(data)/nrow) - if(length(data) == 1 && data == 0 && !identical(sparse,FALSE)) { - if(is.null(sparse)) sparse <- TRUE + if(length(data) == 1 && is0(data) && !identical(sparse, FALSE)) { + ## Matrix(0, ...) : always sparse unless "sparse = FALSE": + if(is.null(sparse)) sparse1 <- sparse <- TRUE ## will be sparse: do NOT construct full matrix! data <- new(if(is.numeric(data)) "dgTMatrix" else if(is.logical(data)) "lgTMatrix" else @@ -135,7 +129,9 @@ dimnames(data) <- dimnames } doDN <- FALSE - } + } else if(!missing(nrow) || !missing(ncol)) + warning("'nrow', 'ncol', etc, are disregarded for matrix 'data'") + ## 'data' is now a "matrix" or "Matrix" if (doDN && !is.null(dimnames)) dimnames(data) <- dimnames @@ -148,11 +144,9 @@ if(isDiag) isDiag <- isDiagonal(data) -### TODO: Compare with as.Matrix() and its tests in ./dgeMatrix.R - ## Find proper matrix class 'cl' cl <- - if(isDiag) + if(isDiag && !isTRUE(sparse1)) "diagonalMatrix" # -> will automatically check for type else { ## consider it's type @@ -193,7 +187,7 @@ function(x, y) callGeneric(x, as.matrix(y))) setMethod("%*%", signature(x = "numeric", y = "Matrix"), - function(x, y) callGeneric(rbind(x), y)) + function(x, y) callGeneric(matrix(x, nrow = 1, byrow=TRUE), y)) setMethod("crossprod", signature(x = "Matrix", y = "numeric"), function(x, y = NULL) callGeneric(x, as.matrix(y))) @@ -207,8 +201,17 @@ setMethod("tcrossprod", signature(x = "numeric", y = "Matrix"), function(x, y = NULL) callGeneric(as.matrix(x), y)) +## maybe not optimal +setMethod("solve", signature(a = "Matrix", b = "missing"), + function(a, b, ...) solve(a, Diagonal(nrow(a)))) + setMethod("solve", signature(a = "Matrix", b = "numeric"), function(a, b, ...) callGeneric(a, as.matrix(b))) +## when no sub-class method is found, bail out +setMethod("solve", signature(a = "Matrix", b = "matrix"), + function(a, b, ...) .bail.out.2("solve", class(a), "matrix")) +setMethod("solve", signature(a = "Matrix", b = "Matrix"), + function(a, b, ...) .bail.out.2("solve", class(a), class(b))) ## bail-out methods in order to get better error messages setMethod("%*%", signature(x = "Matrix", y = "Matrix"), @@ -242,32 +245,20 @@ Y <- as(Y, "matrix") ; Matrix(callGeneric()) }) +## FIXME: All of these should never be called +setMethod("chol", signature(x = "Matrix"), + function(x, pivot = FALSE) .bail.out.1(.Generic, class(x))) +setMethod("determinant", signature(x = "Matrix"), + function(x, logarithm = TRUE) .bail.out.1(.Generic, class(x))) + setMethod("diag", signature(x = "Matrix"), function(x, nrow, ncol) .bail.out.1(.Generic, class(x))) setMethod("t", signature(x = "Matrix"), function(x) .bail.out.1(.Generic, class(x))) ## Group Methods -setMethod("+", signature(e1 = "Matrix", e2 = "missing"), function(e1) e1) -## "fallback": -setMethod("-", signature(e1 = "Matrix", e2 = "missing"), - function(e1) { - warning("inefficient method used for \"- e1\"") - 0-e1 - }) - -## bail-outs: -setMethod("Compare", signature(e1 = "Matrix", e2 = "Matrix"), - function(e1, e2) { - d <- dimCheck(e1,e2) - .bail.out.2(.Generic, class(e1), class(e2)) - }) -setMethod("Compare", signature(e1 = "Matrix", e2 = "ANY"), - function(e1, e2) .bail.out.2(.Generic, class(e1), class(e2))) -setMethod("Compare", signature(e1 = "ANY", e2 = "Matrix"), - function(e1, e2) .bail.out.2(.Generic, class(e1), class(e2))) - +##-> see ./Ops.R ### -------------------------------------------------------------------------- ### @@ -301,30 +292,55 @@ function(x,i,j, drop) stop("invalid or not-yet-implemented 'Matrix' subsetting")) -## "logical *vector* indexing, such as M [ M >= 10 ] : +## logical indexing, such as M[ M >= 7 ] *BUT* also M[ M[,1] >= 3,], +## The following is *both* for M [ ] +## and also for M [ , ] +.M.sub.i.logical <- function (x, i, j, drop) +{ + nA <- nargs() + if(nA == 2) { ## M [ M >= 7 ] + ## FIXME: when both 'x' and 'i' are sparse, this can be very inefficient + as(x, geClass(x))@x[as.vector(i)] + ## -> error when lengths don't match + } else if(nA == 3) { ## M [ M[,1, drop=FALSE] >= 7, ] + stop("not-yet-implemented 'Matrix' subsetting") ## FIXME + + } else stop("nargs() = ", nA, + ". Extraneous illegal arguments inside '[ .. ]' ?") +} setMethod("[", signature(x = "Matrix", i = "lMatrix", j = "missing", drop = "ANY"), - function (x, i, j, drop) { - as(x, geClass(x))@x[as.vector(i)] - ## -> error when lengths don't match - }) - -## FIXME: The following is good for M [ ] -## *BUT* it also triggers for M [ , ] where it is *WRONG* -## using nargs() does not help: it gives '3' for both cases -if(FALSE) + .M.sub.i.logical) setMethod("[", signature(x = "Matrix", i = "logical", j = "missing", drop = "ANY"), - function (x, i, j, drop) { - ## DEBUG - cat("[(Matrix,i,..): nargs=", nargs(),"\n") - as(x, geClass(x))@x[i] }) + .M.sub.i.logical) -## "FIXME:" -## How can we get at A[ ij ] where ij is (i,j) 2-column matrix? -## and A[ LL ] where LL is a logical *vector* -## -> [.data.frame uses nargs() - can we do this in the *generic* ? +## A[ ij ] where ij is (i,j) 2-column matrix : +.M.sub.i.2col <- function (x, i, j, drop) +{ + nA <- nargs() + if(nA == 2) { ## M [ cbind(ii,jj) ] + if(!is.integer(nc <- ncol(i))) + stop("'i' has no integer column number", + " should never happen; please report") + if(is.logical(i)) + return(.M.sub.i.logical(x,i,j,drop)) + else if(!is.numeric(i) || nc != 2) + stop("such indexing must be by logical or 2-column numeric matrix") + m <- nrow(i) + if(m == 0) return(vector(mode = .type.kind[.M.kind(x)])) + ## else + i1 <- i[,1] + i2 <- i[,2] + ## potentially inefficient -- FIXME -- + unlist(lapply(seq_len(m), function(j) x[i1[j], i2[j]])) + + } else stop("nargs() = ", nA, + ". Extraneous illegal arguments inside '[ .. ]' ?") +} +setMethod("[", signature(x = "Matrix", i = "matrix", j = "missing"),# drop="ANY" + .M.sub.i.2col) ### "[<-" : ----------------- @@ -333,49 +349,87 @@ setReplaceMethod("[", signature(x = "Matrix", i = "missing", j = "missing", value = "ANY"),## double/logical/... function (x, value) { - x@x <- value - validObject(x)# check if type and lengths above match - x + ## Fails for 'nMatrix' ... FIXME : make sure have method there + x@x <- rep(value, length = length(x@x)) + validObject(x)# check if type and lengths above match + x }) -## Method for all 'Matrix' kinds (rather than incomprehensible error messages); +## A[ ij ] <- value, where ij is (i,j) 2-column matrix : +## ---------------- The cheap general method --- FIXME: provide special ones +.M.repl.i.2col <- function (x, i, j, value) +{ + nA <- nargs() + if(nA == 3) { ## M [ cbind(ii,jj) ] <- value + if(!is.integer(nc <- ncol(i))) + stop("'i' has no integer column number", + " should never happen; please report") + else if(!is.numeric(i) || nc != 2) + stop("such indexing must be by logical or 2-column numeric matrix") + if(is.logical(i)) { + message(".M.repl.i.2col(): drop 'matrix' case ...") + i <- c(i) # drop "matrix" + return( callNextMethod() ) + } + if(!is.integer(i)) storage.mode(i) <- "integer" + if(any(i < 0)) + stop("negative values are not allowed in a matrix subscript") + if(any(is.na(i))) + stop("NAs are not allowed in subscripted assignments") + if(any(i0 <- (i == 0))) # remove them + i <- i[ - which(i0, arr.ind = TRUE)[,"row"], ] + ## now have integer i >= 1 + m <- nrow(i) + ## mod.x <- .type.kind[.M.kind(x)] + if(length(value) > 0 && m %% length(value) != 0) + warning("number of items to replace is not a multiple of replacement length") + ## recycle: + value <- rep(value, length = m) + i1 <- i[,1] + i2 <- i[,2] + ## inefficient -- FIXME -- (also loses "symmetry" unnecessarily) + for(k in seq_len(m)) + x[i1[k], i2[k]] <- value[k] + + x + } else stop("nargs() = ", nA, + ". Extraneous illegal arguments inside '[ .. ]' ?") +} + +setReplaceMethod("[", signature(x = "Matrix", i = "matrix", j = "missing", + value = "replValue"), + .M.repl.i.2col) + + +setReplaceMethod("[", signature(x = "Matrix", i = "ANY", j = "ANY", + value = "Matrix"), + function (x, i, j, value) { +### *TEMPORARY* diagnostic output: +## cat("[i,j] <- :\n = x :") +## str(x) +## cat(" = value :") +## str(value) +## cat("i :"); if(!missing(i)) str(i) else cat("\n") +## cat("j :"); if(!missing(j)) str(j) else cat("\n") + + callGeneric(x=x, i=i, j=j, value = as.vector(value)) + }) +setReplaceMethod("[", signature(x = "Matrix", i = "ANY", j = "ANY", + value = "Matrix"), + function (x, i, j, value) + callGeneric(x=x, i=i, j=j, value = as.vector(value))) + +setReplaceMethod("[", signature(x = "Matrix", i = "ANY", j = "ANY", + value = "matrix"), + function (x, i, j, value) + callGeneric(x=x, i=i, j=j, value = c(value))) + ## (ANY,ANY,ANY) is used when no `real method' is implemented : setReplaceMethod("[", signature(x = "Matrix", i = "ANY", j = "ANY", value = "ANY"), function (x, i, j, value) { if(!is.atomic(value)) - stop("RHS 'value' must match matrix class ", class(x)) + stop(sprintf("RHS 'value' (class %s) matches 'ANY', but must match matrix class %s", + class(value),class(x))) else stop("not-yet-implemented 'Matrix[<-' method") }) - - -## The trivial methods : -setMethod("cbind2", signature(x = "Matrix", y = "NULL"), - function(x, y) x) -setMethod("cbind2", signature(x = "Matrix", y = "missing"), - function(x, y) x) -setMethod("cbind2", signature(x = "NULL", y="Matrix"), - function(x, y) x) - -setMethod("rbind2", signature(x = "Matrix", y = "NULL"), - function(x, y) x) -setMethod("rbind2", signature(x = "Matrix", y = "missing"), - function(x, y) x) -setMethod("rbind2", signature(x = "NULL", y="Matrix"), - function(x, y) x) - -## Makes sure one gets x decent error message for the unimplemented cases: -setMethod("cbind2", signature(x = "Matrix", y = "Matrix"), - function(x, y) { - rowCheck(x,y) - stop(gettextf("cbind2() method for (%s,%s) not-yet defined", - class(x), class(y))) - }) - -## Use a working fall back {particularly useful for sparse}: -## FIXME: implement rbind2 via "cholmod" for C* and Tsparse ones -setMethod("rbind2", signature(x = "Matrix", y = "Matrix"), - function(x, y) { - colCheck(x,y) - t(cbind2(t(x), t(y))) - })