# R/Auxiliaries.R In Matrix: Sparse and Dense Matrix Classes and Methods

#### Documented in aggregateTanyDuplicatedTasUniqueTcolScale.diag.dsCdiagN2U.diagN2U.diagU2NdimScaledmpermdrop0isUniqueTmat2tripletrowScale.tCRT.validateCsparse

#### "Namespace private" Auxiliaries  such as method functions
#### (called from more than one place --> need to be defined early)

is0  <- function(x) !(is.na(x) | x)
isN0 <- function(x)   is.na(x) | x
is1  <- function(x)  !is.na(x) & x == 1L
isN1 <- function(x)   is.na(x) | x != 1L

allTrue  <- function(x)
!is.na(a <- all( x)) &&  a
allFalse <- function(x)
if(is.atomic(x)) .Call(R_all0, x) else !is.na(a <- any( x)) && !a
all0     <- function(x)
if(is.atomic(x)) .Call(R_all0, x) else !is.na(a <- all(!x)) &&  a
anyFalse <- function(x)
if(is.atomic(x)) .Call(R_any0, x) else !is.na(a <- all( x)) && !a
any0     <- function(x)
if(is.atomic(x)) .Call(R_any0, x) else !is.na(a <- any(!x)) &&  a

## NB: change to using 'typeof' when we define iMatrix
as1 <- function(x, mode. = mode(x))
switch(mode.,
"logical" = TRUE,
"integer" =   1L,
"double"  =     ,
"numeric" =    1,
"complex" = 1+0i,
stop(gettextf("invalid mode \"%s\"", mode.), domain = NA))

## NB: change to using 'typeof' when we define iMatrix
as0 <- function(x, mode. = mode(x))
switch(mode.,
"logical" = FALSE,
"integer" =    0L,
"double"  =      ,
"numeric" =     0,
"complex" =  0+0i,
stop(gettextf("invalid mode \"%s\"", mode.), domain = NA))

.bail.out.2 <- function(name, cl1, cl2)
stop(gettextf("%s(<%s>, <%s>) is not yet implemented; ask maintainer(\"%s\") to implement the missing method",
name, cl1[1L], cl2[1L], "Matrix"),
call. = FALSE, domain = NA)

Matrix.verbose <- function()
getOption("Matrix.verbose", .MatrixEnv[["verbose"]])

Matrix.warn <- function()
getOption("Matrix.warn", .MatrixEnv[["warn"]])

## MJ: or maybe a general one?
if(FALSE) {
Matrix.getOption <- function(x, default = .Matrix.Env[[x]])
getOption(paste0("Matrix.", x), default)
} ## MJ

Matrix.message <- function(..., .M.level = 1, call. = FALSE, domain = NULL) {
if(Matrix.verbose() >= .M.level) {
m <-
if((w <- Matrix.warn()) < 1)
function(..., call., domain) message(..., domain = domain)
else if(w < 2)
warning
else stop
m(..., call. = call., domain = domain)
}
}

## MJ: Implement forceTriangular() and export this and that?
## MJ: Notably this provides a model for (maybe, in the future) allowing
## forceSymmetric(<non-square>) ... by truncating the "too long" dimension.
forceDiagonal <- function(x, diag = NA_character_) {
y <- diag(x, names = FALSE) # FIXME? don't allocate if diag == "U"
cl <- switch(typeof(y),
logical =
{
if(is.na(diag))
diag <-
if(!is.na(a <- all(y        )) && a) "U" else "N"
if(.M.kind(x) == "n") "ndiMatrix" else "ldiMatrix"
},
integer =
{
if(is.na(diag))
diag <-
if(!is.na(a <- all(y == 1L  )) && a) "U" else "N"
storage.mode(y) <- "double"
"ddiMatrix"
},
## integer =
##     {
##         if(is.na(diag))
##             diag <-
##             if(!is.na(a <- all(y == 1L  )) && a) "U" else "N"
##         "idiMatrix"
##     },
double =
{
if(is.na(diag))
diag <-
if(!is.na(a <- all(y == 1   )) && a) "U" else "N"
"ddiMatrix"
},
complex =
stop(gettextf("complex %s not yet implemented",
"diagonalMatrix"),
domain = NA),
## complex =
##     {
##         if(is.na(diag))
##             diag <-
##             if(!is.na(a <- all(y == 1+0i)) && a) "U" else "N"
##         "zdiMatrix"
##     },
stop(gettextf("cannot coerce matrix of type \"%s\" to %s",
typeof(y), "diagonalMatrix"),
domain = NA))
n <- length(y)
d <- dim(x)
dn <- dimnames(x) %||% list(NULL, NULL)
if(any(d > n)) {
d <- c(n, n)
w <- if(d[1L] > n) 1L else 2L
if(!is.null(dnw <- dn[[w]]))
dn[[w]] <- dnw[seq_len(n)]
}
new(cl, Dim = d, Dimnames = dn, diag = diag,
x = if(diag == "N") y else y[FALSE])
}

.tCRT <- function(x, lazy = TRUE) .Call(R_sparse_transpose, x, lazy)

.drop0 <- function(x, tol = 0, isM = TRUE) {
if(isM)
return(.Call(R_sparse_drop0, x, tol))
## TODO: write sparseVector code in C and respecting 'tol'
if(.M.kind(x) == "n")
return(x)
x.x <- x@x
k <- which(is.na(x.x) | x.x)
if(length(k)) {
x@i <- x@i[k]
x@x <- x.x[k]
}
x
}

drop0 <- function(x, tol = 0, is.Csparse = NA, give.Csparse = TRUE) {
tryCoerce <-
if(give.Csparse)
is.na(is.Csparse) || !is.Csparse
else if(.M.kind(x) != "n")
!.isCRT(x)
else if(all(.M.repr(x) != c("C", "R", "T", "i")))
TRUE
else return(x) # n[gst][CRT]Matrix or indMatrix
if(tryCoerce)
x <- if(isS4(x)) .M2C(x) else .m2sparse.checking(x, ".", "C")
.Call(R_sparse_drop0, x, as.double(tol))
}

indDiag <- function(n, upper = TRUE, packed = FALSE)
.Call(R_index_diagonal, n, packed, upper)

indTri <- function(n, upper = TRUE, diag = FALSE, packed = FALSE)
.Call(R_index_triangle, n, packed, upper, diag)

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

##' the workhorse for non0ind(.), but occasionally used directly
non0.i <- function(M, cM = class(M), uniqT = TRUE) {
cld <- getClassDef(cM)
if(extends(cld, "CsparseMatrix"))
.Call(compressed_non_0_ij, M, TRUE)
else if(extends(cld, "RsparseMatrix"))
.Call(compressed_non_0_ij, M, FALSE)
else if(extends(cld, "TsparseMatrix")) {
if(uniqT && !isUniqueT(M))
.Call(compressed_non_0_ij, .M2C(M), TRUE)
else cbind(M@i, M@j, deparse.level = 0L)
} else if(extends(cld, "diagonalMatrix")) {
i <- seq.int(from = 0L, length.out = M@Dim[1L])
if(M@diag == "N")
i <- i[isN0(M@x)]
cbind(i, i, deparse.level = 0L)
} else if(extends(cld, "indMatrix")) {
perm <- M@perm
i <- seq.int(from = 0L, length.out = length(perm))
if(M@margin == 1L)
cbind(i, perm - 1L, deparse.level = 0L)
else cbind(perm - 1L, i, deparse.level = 0L)
} else
stop(gettextf("non0.i() not yet implemented for class %s", dQuote(cM)),
domain = NA)
}

##' the "more versatile / user" function (still not exported):
non0ind <- function(x, cld = getClassDef(class(x)),
uniqT = TRUE, xtendSymm = TRUE, check.Udiag = TRUE)
{
if(is.numeric(x)) {
n <- length(x)
return(if(n == 0L)
integer(0L)
else if(is.matrix(x))
arrayInd(seq_len(n)[isN0(x)], dim(x)) - 1L
else (0:(n-1L))[isN0(x)])
}
stopifnot(extends(cld, "sparseMatrix"))

ij <- non0.i(x, cld, uniqT = uniqT)
if(xtendSymm && extends(cld, "symmetricMatrix")) {
## also get "other" triangle, but not the diagonal again
notdiag <- ij[, 1L] != ij[, 2L]
rbind(ij, ij[notdiag, 2:1], deparse.level = 0L)
} else if(check.Udiag && extends(cld, "triangularMatrix") &&
x@diag == "U") {
i <- seq.int(from = 0L, length.out = x@Dim[1L])
rbind(ij, cbind(i, i, deparse.level = 0L), deparse.level = 0L)
} else ij
}

##' Decode "encoded" (i,j) indices back to  cbind(i,j)
##' This is the inverse of encodeInd(.)
##'
##' @title Decode "Encoded" (i,j) Indices
##' @param code integer in 0:((n x m - 1)  <==> encodeInd(.) result
##' @param nr the number of rows
##' @return
##' @author Martin Maechler
decodeInd <- function(code, nr)
cbind(as.integer(code %% nr), as.integer(code %/% nr),
deparse.level = 0L)

complementInd <- function(ij, dim, orig1=FALSE, checkBnds=FALSE) {
## Purpose: Compute the complement of the 2-column 0-based ij-matrix
##		but as 1-based indices
n <- prod(dim)
if(n == 0L)
return(integer(0L))
seq_len(n)[-(1L + .Call(m_encodeInd, ij, dim, orig1, checkBnds))]
}

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

anyDuplicatedT <- function(x, ...) {
mn <- prod(d <- x@Dim)
if(mn <= .Machine\$integer.max)
anyDuplicated.default(          x@j  * d[1L] + x@i, ...)
else if(mn <= 0x1p+53)
anyDuplicated.default(as.double(x@j) * d[1L] + x@i, ...)
else anyDuplicated.default(.mapply(c, list(x@i, x@j), NULL), ...)
}

isUniqueT <- function(x, byrow = FALSE, isT = is(x, "TsparseMatrix"))
isT && !is.unsorted(if(byrow) order(x@i, x@j) else order(x@j, x@i),
strictly = TRUE)

asUniqueT <- function(x, byrow = FALSE, isT = is(x, "TsparseMatrix"))
if(isUniqueT(x, byrow, isT)) x else .M2T(if(byrow) .M2R(x) else .M2C(x))

aggregateT <- function(x) .Call(Tsparse_aggregate, x)

mat2triplet <- function(x, uniqT = FALSE) {
T <- as(x, "TsparseMatrix")
if(uniqT)
T <- asUniqueT(T, isT = TRUE)
if(is(T, "nsparseMatrix"))
list(i = T@i + 1L, j = T@j + 1L)
else list(i = T@i + 1L, j = T@j + 1L, x = T@x)
}

.validateCsparse <- function(x, sort.if.needed = FALSE) {
if(sort.if.needed)
.Call(CsparseMatrix_validate_maybe_sorting, x)
else .Call(CsparseMatrix_validate, x)
}

dmperm <- function(x, nAns = 6L, seed = 0L) {
## NB: result is determined entirely by nonzero pattern of 'x'
stopifnot(length(nAns <- as.integer(nAns)) == 1L, any(nAns == 2L * (1L:3L)),
length(seed <- as.integer(seed)) == 1L, any(seed == -1L:1L))
x <- if(isS4(x)) .M2gen(.M2C(x), "n") else .m2sparse(x, "ngC")
.Call(Csparse_dmperm, x, nAns, seed)
}

## (matrix|denseMatrix)->denseMatrix as similar as possible to "target"
as_denseClass <- function(x, cl, cld = getClassDef(cl)) {
kind <- .M.kind(x)
cl <- .M.nonvirtual(new(cld))
if(cl == "indMatrix")
cl <- "ngeMatrix"
symmetric <- substr(cl, 2L, 2L) == "s" && isSymmetric(x)
triangular <- !symmetric &&
substr(cl, 2L, 2L) == "t" && (it <- isTriangular(x))
packed <- (symmetric || triangular) && substr(cl, 3L, 3L) == "p"
if(isS4(x)) {
r <- if(symmetric || triangular)
.M2kind(if(symmetric)
forceSymmetric(x)
else if(attr(it, "kind") == "U")
triu(x)
else tril(x),
kind)
else .M2gen(x, kind)
if(packed) pack(r) else r
}
else if(symmetric)
.m2dense(x, paste0(kind, "s", if(packed) "p" else "y"), "U", NULL)
else if(triangular)
.m2dense(x, paste0(kind, "t", if(packed) "p" else "r"), attr(it, "kind"), "N")
else .m2dense(x, paste0(kind, "ge"))
}

## (matrix|sparseMatrix)->CsparseMatrix as similar as possible to "target"
as_CspClass <- function(x, cl, cld = getClassDef(cl)) {
kind <- .M.kind(x)
cl <- .M.nonvirtual(new(cld))
if(cl == "indMatrix")
cl <- "ngeMatrix"
symmetric <- substr(cl, 2L, 2L) == "s" && isSymmetric(x)
triangular <- !symmetric &&
substr(cl, 2L, 2L) == "t" && (it <- isTriangular(x))
if(isS4(x)) {
r <- if(symmetric || triangular)
.M2kind(if(symmetric)
forceSymmetric(x)
else if(attr(it, "kind") == "U")
triu(x)
else tril(x),
kind)
else .M2gen(x, kind)
.M2C(r)
}
else if(symmetric)
.m2sparse(x, paste0(kind, "sC"), "U", NULL)
else if(triangular)
.m2sparse(x, paste0(kind, "tC"), attr(it, "kind"), "N")
else .m2sparse(x, paste0(kind, "gC"))
}

## as(<[Mm]atrix>, <non-unit triangular CsparseMatrix>)
asCspN <- function(x) .Call(R_sparse_diag_U2N, as(x, "CsparseMatrix"))

diagU2N <- function (x, cl = getClassDef(class(x)), checkDense = FALSE) {
if(extends(cl, "triangularMatrix") && x@diag == "U")
.diagU2N(x, cl = cl, checkDense = checkDense)
else x
}

.diagU2N <- function(x, cl = getClassDef(class(x)), checkDense = FALSE) {
if(!checkDense && extends(cl, "denseMatrix"))
x <- as(x, "CsparseMatrix")
..diagU2N(x)
}

..diagU2N <- function(x) {
diag(x) <- TRUE
x
}

diagN2U <- function(x, cl = getClassDef(class(x)), checkDense = FALSE) {
if(extends(cl, "triangularMatrix") && x@diag == "N")
.diagN2U(x, cl = cl, checkDense = checkDense)
else x
}

.diagN2U <- function(x, cl = getClassDef(class(x)), checkDense = FALSE) {
if(!checkDense & (isDense <- extends(cl, "denseMatrix")))
..diagN2U(as(x, "CsparseMatrix"), sparse = TRUE)
else ..diagN2U(x, sparse = !isDense)
}

..diagN2U <- function(x, sparse) {
if(sparse && x@Dim[1L] > 0L)
x <- switch(x@uplo,
U = .Call(R_sparse_band, x, 1L, NULL),
L = .Call(R_sparse_band, x, NULL, -1L))
x@diag <- "U"
x
}

## Caches 'value' in the 'factors' slot of 'x' { NOT a copy of 'x' ! }
## and returns 'value'
.set.factor <- function(x, name, value, warn.no.slot = FALSE)
.Call(R_set_factor, x, name, value, warn.no.slot)

##' Compute the three "parts" of two sets:
.setparts <- function(x, y) {
n1 <- length(m1 <- match(x, y, 0L))
n2 <- length(m2 <- match(y, x, 0L))
ix <- which(m1 == 0L)
iy <- which(m2 == 0L)
list(x.only = x[ix], ix.only = ix, mx = m1,
y.only = y[iy], iy.only = iy, my = m2,
int = if(n1 < n2) y[m1] else x[m2])
}

##' *Only* to be used as function in
##'    setMethod.("Compare", ...., .Cmp.swap)  -->  ./Ops.R  & ./diagMatrix.R
.Cmp.swap <- function(e1, e2) {
## "swap RHS and LHS" and use the method below:
switch(.Generic,
"==" =,
"!=" = callGeneric(e2, e1),
"<"	= e2 >	e1,
"<=" = e2 >= e1,
">"	= e2 <	e1,
">=" = e2 <= e1)
}

.diag.dsC <- function(x, Chx = Cholesky(x, LDL = TRUE), res.kind = "diag") {
if(!missing(Chx))
stopifnot(.CHM.is.LDL(Chx), is.integer(Chx@p), is.double(Chx@x))
.Call(dtCMatrix_diag, Chx, res.kind)
##    ^^^^^^^^^^^^^^ in ../src/Csparse.c
}

dimScale <- function(x, d1 = sqrt(1/diag(x, names = FALSE)), d2 = d1) {
dim.x <- dim(x)
D1 <- Diagonal(n = dim.x[1L], x = d1)
D2 <- if(missing(d2)) D1 else Diagonal(n = dim.x[2L], x = d2)
y <- D1 %*% x %*% D2 # inefficient for symmetricMatrix 'x', but "general"
if(isS4(x) && is(x, "symmetricMatrix") && identical(d1, d2))
y <- forceSymmetric(y, x@uplo)
if(is.list(dn <- dimnames(x)))
y@Dimnames <- dn
y
}

rowScale <- function(x, d) {
y <- Diagonal(n = nrow(x), x = d) %*% x
if(is.list(dn <- dimnames(x)))
y@Dimnames <- dn
y
}

colScale <- function(x, d) {
y <- x %*% Diagonal(n = ncol(x), x = d)
if(is.list(dn <- dimnames(x)))
y@Dimnames <- dn
y
}


## Try the Matrix package in your browser

Any scripts or data that you put into this service are public.

Matrix documentation built on May 29, 2024, 1:20 a.m.