SCM Repository

[matrix] View of /pkg/R/Auxiliaries.R
 [matrix] / pkg / R / Auxiliaries.R

View of /pkg/R/Auxiliaries.R

Thu Sep 28 15:31:17 2006 UTC (13 years, 5 months ago) by maechler
File size: 19724 byte(s)
"Compare" for "dMatrix"; plus a few small ones
#### "Namespace private" Auxiliaries  such as method functions
#### (called from more than one place --> need to be defined early)

.isR_24 <- (paste(R.version$major, R.version$minor, sep=".") >= "2.4")

## Need to consider NAs ;  "== 0" even works for logical & complex:
is0  <- function(x) !is.na(x) & x == 0
isN0 <- function(x)  is.na(x) | x != 0
all0 <- function(x) !any(is.na(x)) && all(x == 0)

allTrue  <- function(x) !any(is.na(x)) && all(x)
allFalse <- function(x) !any(is.na(x)) && !any(x)

## For %*% (M = Matrix; v = vector (double or integer {complex maybe?}):
.M.v <- function(x, y) callGeneric(x, as.matrix(y))
.v.M <- function(x, y) callGeneric(rbind(x), y)

.M.DN <- function(x) if(!is.null(dn <- dimnames(x))) dn else list(NULL,NULL)

.has.DN <- ## has non-trivial Dimnames slot?
function(x) !identical(list(NULL,NULL), x@Dimnames)

.bail.out.1 <- function(fun, cl) {
stop(gettextf('not-yet-implemented method for %s(<%s>)', fun, cl),
call. = FALSE)
}
.bail.out.2 <- function(fun, cl1, cl2) {
stop(gettextf('not-yet-implemented method for %s(<%s>, <%s>)',
fun, cl1, cl2), call. = FALSE)
}

## This should be done in C and be exported by 'methods':  [FIXME - ask JMC ]
copyClass <- function(x, newCl, sNames =
intersect(slotNames(newCl), slotNames(x))) {
r <- new(newCl)
for(n in sNames)
slot(r, n) <- slot(x, n)
r
}

## chol() via "dpoMatrix"
cholMat <- function(x, pivot, ...) {
px <- as(x, "dpoMatrix")
if (isTRUE(validObject(px, test=TRUE))) chol(px)
else stop("'x' is not positive definite -- chol() undefined.")
}

dimCheck <- function(a, b) {
da <- dim(a)
db <- dim(b)
if(any(da != db))
stop(gettextf("Matrices must have same dimensions in %s",
deparse(sys.call(sys.parent()))),
call. = FALSE)
da
}

dimNamesCheck <- function(a, b) {
## assume dimCheck() has happened before
nullDN <- list(NULL,NULL)
h.a <- !identical(nullDN, dna <- dimnames(a))
h.b <- !identical(nullDN, dnb <- dimnames(b))
if(h.a || h.b) {
if (!h.b) dna
else if(!h.a) dnb
else { ## both have non-trivial dimnames
r <- dna # "default" result
for(j in 1:2) {
dn <- dnb[[j]]
if(is.null(r[[j]]))
r[[j]] <- dn
else if (!is.null(dn) && any(r[[j]] != dn))
warning(gettextf("dimnames [%d] mismatch in %s", j,
deparse(sys.call(sys.parent()))),
call. = FALSE)
}
r
}
}
else
nullDN
}

rowCheck <- function(a, b) {
da <- dim(a)
db <- dim(b)
if(da[1] != db[1])
stop(gettextf("Matrices must have same number of rows in %s",
deparse(sys.call(sys.parent()))),
call. = FALSE)
## return the common nrow()
da[1]
}

colCheck <- function(a, b) {
da <- dim(a)
db <- dim(b)
if(da[2] != db[2])
stop(gettextf("Matrices must have same number of columns in %s",
deparse(sys.call(sys.parent()))),
call. = FALSE)
## return the common ncol()
da[2]
}

## Note: !isPacked(.)  i.e. full' still contains
## ----  "*sy" and "*tr" which have "undefined" lower or upper part
isPacked <- function(x)
{
## Is 'x' a packed (dense) matrix ?
is(x, "denseMatrix") &&
any("x" == slotNames(x)) && length(x@x) < prod(dim(x))
}

emptyColnames <- function(x)
{
## Useful for compact printing of (parts) of sparse matrices
## possibly	 dimnames(x) "==" NULL :
dimnames(x) <- list(dimnames(x)[[1]], rep("", dim(x)[2]))
x
}

### TODO:  write in C and port to base (or 'utils') R
indTri <- function(n, upper = TRUE) {
## == which(upper.tri(diag(n)) or
##	  which(lower.tri(diag(n)) -- but much more efficiently for largish 'n'
stopifnot(length(n) == 1, n == (n. <- as.integer(n)), (n <- n.) >= 0)
if(n <= 2)
return(if(n == 2) as.integer(if(upper) n+1 else n) else integer(0))
## First, compute the 'diff(.)'  fast.  Use integers
one <- 1:1 ; two <- 2:2
n1 <- n - one
n2 <- n1 - one
r <- rep.int(one, n*n1/two - one)
r[cumsum(if(upper) 1:n2 else c(n1, if(n >= 4) n2:two))] <- if(upper) n:3 else 3:n
## now have "dliu" difference; revert to "liu":
cumsum(c(if(upper) n+one else two, r))
}

prTriang <- function(x, digits = getOption("digits"),
maxp = getOption("max.print"),
justify = "none", right = TRUE)
{
## modeled along stats:::print.dist
upper <- x@uplo == "U"

m <- as(x, "matrix")
cf <- format(m, digits = digits, justify = justify)
if(upper)
cf[row(cf) > col(cf)] <- "."
else
cf[row(cf) < col(cf)] <- "."
if(.isR_24)
print(cf, quote = FALSE, right = right, max = maxp)
else print(cf, quote = FALSE, right = right)
invisible(x)
}

prMatrix <- function(x, digits = getOption("digits"),
maxp = getOption("max.print")) {
d <- dim(x)
cl <- class(x)
cat(sprintf('%d x %d Matrix of class "%s"\n', d[1], d[2], cl))
if(prod(d) <= maxp) {
if(is(x, "triangularMatrix"))
prTriang(x, digits = digits, maxp = maxp)
else {
if(.isR_24)
print(as(x, "matrix"), digits = digits, max = maxp)
else print(as(x, "matrix"), digits = digits)
}
}
else { ## d[1] > maxp / d[2] >= nr :
m <- as(x, "matrix")
nr <- maxp %/% d[2]
n2 <- ceiling(nr / 2)
cat("\n ..........\n\n")
print(tail(m, max(1, nr - n2)))
}
## DEBUG: cat("str(.):\n") ; str(x)
invisible(x)# as print() S3 methods do
}

nonFALSE <- function(x) {
## typically used for lMatrices:  (TRUE,NA,FALSE) |-> (TRUE,FALSE)
if(any(ix <- is.na(x))) x[ix] <- TRUE
x
}

nz.NA <- function(x, na.value) {
## Non-Zeros of x
## na.value: TRUE: NA's give TRUE, they are not 0
##             NA: NA's are not known ==> result := NA
##          FALSE: NA's give FALSE, could be 0
stopifnot(is.logical(na.value) && length(na.value) == 1)
if(is.na(na.value)) x != 0
else  if(na.value)	isN0(x)
else		x != 0 & !is.na(x)
}

## Number of non-zeros :
## FIXME? -- make this into a generic function (?)
nnzero <- function(x, na.counted = NA) {
## na.counted: TRUE: NA's are counted, they are not 0
##               NA: NA's are not known (0 or not) ==>  result := NA
##            FALSE: NA's are omitted before counting
cl <- class(x)
if(!extends(cl, "Matrix"))
sum(nz.NA(x, na.counted))
else if(extends(cl, "sparseMatrix"))
## NOTA BENE: The number of *structural* non-zeros {could have other '0'}!
switch(.sp.class(cl),
"CsparseMatrix" = length(x@i),
"TsparseMatrix" = length(x@i),
"RsparseMatrix" = length(x@j))
else ## denseMatrix
sum(nz.NA(as_geClass(x, cl)@x, na.counted))
}

## For sparseness handling
## return a 2-column (i,j) matrix of
## 0-based indices of non-zero entries  :
non0ind <- function(x) {

if(is.numeric(x))
return(if((n <- length(x))) (0:(n-1))[isN0(x)] else integer(0))
## else
stopifnot(is(x, "sparseMatrix"))
non0.i <- function(M) {
if(is(M, "TsparseMatrix"))
return(unique(cbind(M@i,M@j)))
if(is(M, "pMatrix"))
return(cbind(seq(length=nrow(M)), M@perm) - 1:1)
## else:
isC <- any("i" == slotNames(M)) # is Csparse (not Rsparse)
.Call(compressed_non_0_ij, M, isC)
}

if(is(x, "symmetricMatrix")) { # also get "other" triangle
ij <- non0.i(x)
notdiag <- ij[,1] != ij[,2]# but not the diagonals again
rbind(ij, ij[notdiag, 2:1])
}
else if(is(x, "triangularMatrix")) { # check for "U" diag
if(x@diag == "U") {
i <- seq(length = dim(x)[1]) - 1:1
rbind(non0.i(x), cbind(i,i))
} else non0.i(x)
}
else
non0.i(x)
}

## nr= nrow: since  i in {0,1,.., nrow-1}  these are 1:1 "decimal" encodings:
## Further, these map to and from the usual "Fortran-indexing" (but 0-based)
encodeInd <- function(ij, nr) ij[,1] + ij[,2] * nr
decodeInd <- function(code, nr) cbind(code %% nr, code %/% nr)

complementInd <- function(ij, dim)
{
## Purpose: Compute the complement of the 2-column 0-based ij-matrix
##		but as 1-based indices
n <- prod(dim)
if(n == 0) return(integer(0))
ii <- 1:n
ii[-(1 + encodeInd(ij, nr = dim[1]))]
}

unionInd <- function(ij1, ij2) unique(rbind(ij1, ij2))

intersectInd <- function(ij1, ij2, nrow) {
## from 2-column (i,j) matrices where i in {0,.., nrow-1},
## return only the *common* entries
decodeInd(intersect(encodeInd(ij1, nrow),
encodeInd(ij2, nrow)), nrow)
}

WhichintersectInd <- function(ij1, ij2, nrow) {
## from 2-column (i,j) matrices where i \in {0,.., nrow-1},
## find *where*  common entries are in ij1 & ij2
m1 <- match(encodeInd(ij1, nrow), encodeInd(ij2, nrow))
ni <- !is.na(m1)
list(which(ni), m1[ni])
}

### There is a test on this in ../tests/dgTMatrix.R !

uniqTsparse <- function(x, class.x = c(class(x))) {
## Purpose: produce a *unique* triplet representation:
##		by having (i,j) sorted and unique
## -----------------------------------------------------------
## The following is not quite efficient {but easy to program,
## and as() are based on C code  (all of them?)
##
## FIXME: Do it fast for the case where 'x' is already 'uniq'

switch(class.x,
"dgTMatrix" = as(as(x, "dgCMatrix"), "dgTMatrix"),
"dsTMatrix" = as(as(x, "dsCMatrix"), "dsTMatrix"),
"dtTMatrix" = as(as(x, "dtCMatrix"), "dtTMatrix"),
## do we need this for "logical" ones, there's no sum() there!
"lgTMatrix" = as(as(x, "lgCMatrix"), "lgTMatrix"),
"lsTMatrix" = as(as(x, "lsCMatrix"), "lsTMatrix"),
"ltTMatrix" = as(as(x, "ltCMatrix"), "ltTMatrix"),
## do we need this for "logical" ones, there's no sum() there!
"ngTMatrix" = as(as(x, "ngCMatrix"), "ngTMatrix"),
"nsTMatrix" = as(as(x, "nsCMatrix"), "nsTMatrix"),
"ntTMatrix" = as(as(x, "ntCMatrix"), "ntTMatrix"),
## otherwise:
stop("not yet implemented for class ", class.x))
}

## Note: maybe, using
## ----    xj <- .Call(Matrix_expand_pointers, x@p)
## would be slightly more efficient than as( <dgC> , "dgTMatrix")
## but really efficient would be to use only one .Call(.) for uniq(.) !

uniq <- function(x) {
if(is(x, "TsparseMatrix")) uniqTsparse(x) else x
## else:  not 'Tsparse', i.e. "uniquely" represented in any case
}

asTuniq <- function(x) {
if(is(x, "TsparseMatrix")) uniqTsparse(x) else as(x,"TsparseMatrix")
}

if(FALSE) ## try an "efficient" version
uniq_gT <- function(x)
{
## Purpose: produce a *unique* triplet representation:
##		by having (i,j) sorted and unique
## ------------------------------------------------------------------
## Arguments: a "gT" Matrix
stopifnot(is(x, "gTMatrix"))
if((n <- length(x@i)) == 0) return(x)
ii <- order(x@i, x@j)
if(any(ii != 1:n)) {
x@i <- x@i[ii]
x@j <- x@j[ii]
x@x <- x@x[ii]
}
ij <- x@i + nrow(x) * x@j
if(any(dup <- duplicated(ij))) {

}
### We should use a .Call() based utility for this!

}

t_geMatrix <- function(x) {
x@x <- as.vector(t(array(x@x, dim = x@Dim))) # no dimnames here
x@Dim <- x@Dim[2:1]
x@Dimnames <- x@Dimnames[2:1]
## FIXME: how to set factors?
x
}

## t( [dl]trMatrix ) and  t( [dl]syMatrix ) :
t_trMatrix <- function(x) {
x@x <- as.vector(t(as(x, "matrix")))
x@Dim <- x@Dim[2:1]
x@Dimnames <- x@Dimnames[2:1]
x@uplo <- if (x@uplo == "U") "L" else "U"
# and keep x@diag
x
}

fixupDense <- function(m, from) {
if(is(m, "triangularMatrix")) {
m@uplo <- from@uplo
m@diag <- from@diag
} else if(is(m, "symmetricMatrix")) {
m@uplo <- from@uplo
}
m
}

## -> ./ldenseMatrix.R :
l2d_Matrix <- function(from) {
stopifnot(is(from, "lMatrix"))
fixupDense(new(sub("^l", "d", class(from)),
x = as.double(from@x),
Dim = from@Dim, Dimnames = from@Dimnames),
from)
## FIXME: treat 'factors' smartly {not for triangular!}
}

## -> ./ndenseMatrix.R :
n2d_Matrix <- function(from) {
stopifnot(is(from, "nMatrix"))
fixupDense(new(sub("^n", "d", class(from)),
x = as.double(from@x),
Dim = from@Dim, Dimnames = from@Dimnames),
from)
## FIXME: treat 'factors' smartly {not for triangular!}
}
n2l_spMatrix <- function(from) {
stopifnot(is(from, "nMatrix"))
new(sub("^n", "l", class(from)),
##x = as.double(from@x),
Dim = from@Dim, Dimnames = from@Dimnames)
}

if(FALSE)# unused
l2d_meth <- function(x) {
cl <- class(x)
as(callGeneric(as(x, sub("^l", "d", cl))), cl)
}

## return "d" or "l" or "n" or "z"
.M.kind <- function(x, clx = class(x)) {
if(is.matrix(x)) { ## 'old style matrix'
if     (is.numeric(x)) "d"
else if(is.logical(x)) "l" ## FIXME ? "n" if no NA ??
else if(is.complex(x)) "z"
else stop("not yet implemented for matrix w/ typeof ", typeof(x))
}
else if(extends(clx, "dMatrix")) "d"
else if(extends(clx, "nMatrix")) "n"
else if(extends(clx, "lMatrix")) "l"
else if(extends(clx, "zMatrix")) "z"
else if(extends(clx, "pMatrix")) "n" # permutation -> pattern
else stop(" not yet be implemented for ", clx)
}

.M.shape <- function(x, clx = class(x)) {
if(is.matrix(x)) { ## 'old style matrix'
if     (isDiagonal  (x)) "d"
else if(isTriangular(x)) "t"
else if(isSymmetric (x)) "s"
else "g" # general
}
else if(extends(clx, "diagonalMatrix"))  "d"
else if(extends(clx, "triangularMatrix"))"t"
else if(extends(clx, "symmetricMatrix")) "s"
else "g"
}

class2 <- function(cl, kind = "l", do.sub = TRUE) {
## Find "corresponding" class; since pos.def. matrices have no pendant:
if	   (cl == "dpoMatrix") paste(kind, "syMatrix", sep='')
else if(cl == "dppMatrix") paste(kind, "spMatrix", sep='')
else if(do.sub) sub("^d", kind, cl)
else cl
}

geClass <- function(x) {
if     (is(x, "dMatrix")) "dgeMatrix"
else if(is(x, "lMatrix")) "lgeMatrix"
else if(is(x, "nMatrix")) "ngeMatrix"
else if(is(x, "zMatrix")) "zgeMatrix"
else stop("general Matrix class not yet implemented for ",
class(x))
}

.dense.prefixes <- c("d" = "di",
"t" = "tr",
"s" = "sy",
"g" = "ge")

.sparse.prefixes <- c("d" = "t", ## map diagonal to triangular
"t" = "t",
"s" = "s",
"g" = "g")

## Used, e.g. after subsetting: Try to use specific class -- if feasible :
as_dense <- function(x) {
as(x, paste(.M.kind(x), .dense.prefixes[.M.shape(x)], "Matrix", sep=''))
}

.sp.class <- function(x) { ## find and return the "sparseness class"
if(!is.character(x)) x <- class(x)
for(cl in paste(c("C","T","R"), "sparseMatrix", sep=''))
if(extends(x, cl))
return(cl)
## else (should rarely happen)
as.character(NA)
}

as_Csparse <- function(x) {
as(x, paste(.M.kind(x), .sparse.prefixes[.M.shape(x)], "CMatrix", sep=''))
}
as_Rsparse <- function(x) {
as(x, paste(.M.kind(x), .sparse.prefixes[.M.shape(x)], "RMatrix", sep=''))
}
as_Tsparse <- function(x) {
as(x, paste(.M.kind(x), .sparse.prefixes[.M.shape(x)], "TMatrix", sep=''))
}

as_geClass <- function(x, cl) {
if	   (extends(cl, "diagonalMatrix")  && isDiagonal(x))
as(x, cl)
else if(extends(cl, "symmetricMatrix") &&  isSymmetric(x)) {
kind <- .M.kind(x)
as(x, class2(cl, kind, do.sub= kind != "d"))
} else if(extends(cl, "triangularMatrix") && isTriangular(x))
as(x, cl)
else
as(x, paste(.M.kind(x), "geMatrix", sep=''))
}

as_CspClass <- function(x, cl) {
if ((extends(cl, "diagonalMatrix")	&& isDiagonal(x)) ||
(extends(cl, "symmetricMatrix") &&  isSymmetric(x)) ||
(extends(cl, "triangularMatrix")&& isTriangular(x)))
as(x, cl)
else as(x, paste(.M.kind(x), "gCMatrix", sep=''))
}

## -> ./ddenseMatrix.R :
d2l_Matrix <- function(from) {
stopifnot(is(from, "dMatrix"))
fixupDense(new(sub("^d", "l", class(from)), # no need for dClass2 here
Dim = from@Dim, Dimnames = from@Dimnames),
from)
## FIXME: treat 'factors' smartly {not for triangular!}
}

try_as <- function(x, classes, tryAnyway = FALSE) {
if(!tryAnyway && !is(x, "Matrix"))
return(x)
## else
ok <- canCoerce(x, classes[1])
while(!ok && length(classes <- classes[-1])) {
ok <- canCoerce(x, classes[1])
}
if(ok) as(x, classes[1]) else x
}

## For *dense* matrices
isTriMat <- function(object, upper = NA) {
## pretest: is it square?
d <- dim(object)
if(d[1] != d[2]) return(FALSE)
## else slower test
if(!is.matrix(object))
object <- as(object,"matrix")
if(is.na(upper)) {
if(all0(object[lower.tri(object)]))
structure(TRUE, kind = "U")
else if(all0(object[upper.tri(object)]))
structure(TRUE, kind = "L")
else FALSE
} else if(upper)
all0(object[lower.tri(object)])
else ## upper is FALSE
all0(object[upper.tri(object)])
}

## For Csparse matrices
isTriC <- function(x, upper = NA) {
## pretest: is it square?
d <- dim(x)
if(d[1] != d[2]) return(FALSE)
## else
if(d[1] == 0) return(TRUE)
ni <- 1:d[2]
## the row indices split according to column:
ilist <- split(x@i, factor(rep.int(ni, diff(x@p)), levels= ni))
lil <- unlist(lapply(ilist, length), use.names = FALSE)
if(any(lil == 0)) {
pos <- lil > 0
if(!any(pos)) ## matrix of all 0's
return(TRUE)
ilist <- ilist[pos]
ni <- ni[pos]
}
ni0 <- ni - 1:1 # '0-based ni'
if(is.na(upper)) {
if(all(sapply(ilist, max, USE.NAMES = FALSE) <= ni0))
structure(TRUE, kind = "U")
else if(all(sapply(ilist, min, USE.NAMES = FALSE) >= ni0))
structure(TRUE, kind = "L")
else FALSE
} else if(upper) {
all(sapply(ilist, max, USE.NAMES = FALSE) <= ni0)
} else { ## 'lower'
all(sapply(ilist, min, USE.NAMES = FALSE) >= ni0)
}
}

.is.diagonal <- function(object) {
## "matrix" or "denseMatrix" (but not "diagonalMatrix")
d <- dim(object)
if(d[1] != (n <- d[2])) FALSE
else if(is.matrix(object))
## requires that "vector-indexing" works for 'object' :
all0(object[rep(c(FALSE, rep.int(TRUE,n)), length = n^2)])
else ## "denseMatrix" -- packed or unpacked
if(is(object, "generalMatrix")) # "dge", "lge", ...
all0(object@x[rep(c(FALSE, rep.int(TRUE,n)), length = n^2)])
else { ## "dense" but not {diag, general}, i.e. triangular or symmetric:
## -> has 'uplo'  differentiate between packed and unpacked

### .......... FIXME ...............

packed <- isPacked(object)
if(object@uplo == "U") {
} else { ## uplo == "L"
}

### very cheap workaround
all0(as.matrix(object)[rep(c(FALSE, rep.int(TRUE,n)), length = n^2)])
}
}

## FIXME? -- this should also work for "ltT", "ntT", ... :
diagU2N <- function(x)
{
## Purpose: Transform a *unit diagonal* sparse triangular matrix
##	into one with explicit diagonal entries '1'
xT <- as(x, "dgTMatrix")
## leave it as  T* - the caller can always coerce to C* if needed:
new("dtTMatrix", x = xT@x, i = xT@i, j = xT@j, Dim = x@Dim,
Dimnames = x@Dimnames, uplo = x@uplo, diag = "N")
}

## FIXME: this should probably be dropped / replaced by as_Csparse
.as.dgC.Fun <- function(x, na.rm = FALSE, dims = 1) {
x <- as(x, "dgCMatrix")
callGeneric()
}

.as.dgT.Fun <- function(x, na.rm = FALSE, dims = 1) {
## used e.g. inside colSums() etc methods
x <- as(x, "dgTMatrix")
callGeneric()
}

### Fast much simplified version of tapply()
tapply1 <- function (X, INDEX, FUN = NULL, ..., simplify = TRUE) {
sapply(split(X, INDEX), FUN, ..., simplify = simplify, USE.NAMES = FALSE)
}

## tapply.x <- function (X, n, INDEX, FUN = NULL, ..., simplify = TRUE) {
##     tapply1(X, factor(INDEX, 0:(n-1)), FUN = FUN, ..., simplify = simplify)
## }

`