The R Project SVN R-packages

#### Will be sourced by several R scripts in ../tests/

### ------- Part I --  unrelated to "Matrix" classes ---------------

paste0 <- function(...) paste(..., sep = '')

identical3 <- function(x,y,z)   identical(x,y) && identical (y,z)
identical4 <- function(a,b,c,d) identical(a,b) && identical3(b,c,d)

## Make sure errors are signaled
assertError <- function(expr) {
    d.expr <- deparse(substitute(expr))
    t.res <- try(expr, silent = TRUE)
    if(!inherits(t.res, "try-error"))
    stop(d.expr, "\n\t did not give an error", call. = FALSE)
    invisible(t.res)
}

isValid <- function(x, class) validObject(x, test=TRUE) && is(x, class)

is.all.equal3 <- function(x,y,z, tol = .Machine$double.eps^0.5)
    isTRUE(all.equal(x,y, tol=tol)) && isTRUE(all.equal(y,z, tol=tol))

is.all.equal4 <- function(x,y,z,u, tol = .Machine$double.eps^0.5)
    is.all.equal3(x,y,z, tol=tol) && isTRUE(all.equal(z,u, tol=tol))

## A version of all.equal() for the slots
all.slot.equal <- function(x,y, ...) {
    slts <- slotNames(x)
    for(sl in slts) {
        aeq <- all.equal(slot(x,sl), slot(y,sl), ...)
        if(!identical(TRUE, aeq))
            return(paste("slot '",sl,"': ", aeq, sep=''))
    }
    TRUE
}

## The relative error typically returned by all.equal:
relErr <- function(target, current) { ## make this work for 'Matrix'
    ## ==> no mean() ..
    n <- length(current)
    if(length(target) < n)
        target <- rep(target, length.out = n)
    sum(abs(target - current)) / sum(abs(target))
}

## is.R22 <- (paste(R.version$major, R.version$minor, sep=".") >= "2.2")

pkgRversion <- function(pkgname)
    substring(packageDescription(pkgname)[["Built"]], 3,5)


### ------- Part II  -- related to matrices, but *not* "Matrix" -----------

add.simpleDimnames <- function(m) {
    stopifnot(length(d <- dim(m)) == 2)
    dimnames(m) <- list(paste0("r", seq_len(d[1])),
                        paste0("c", seq_len(d[2])))
    m
}

as.mat <- function(m) {
    ## as(., "matrix")  but with no extraneous empty dimnames
    m <- as(m, "matrix")
    if(identical(dimnames(m), list(NULL,NULL)))
    dimnames(m) <- NULL
    m
}

assert.EQ.mat <- function(M, m, tol = if(show) 0 else 1e-15, show=FALSE) {
    ## Purpose: check equality of  'Matrix' M with  'matrix' m
    ## ----------------------------------------------------------------------
    ## Arguments: M: is(., "Matrix") typically {but just needs working as(., "matrix")}
    ##            m: is(., "matrix")
    ##            show: if TRUE, return (and hence typically print) all.equal(...)
    MM <- as.mat(M)                     # as(M, "matrix")
    if(is.logical(MM) && is.numeric(m))
    storage.mode(MM) <- "integer"
    attr(MM, "dimnames") <- attr(m, "dimnames") <- NULL
    if(show) all.equal(MM, m, tol = tol)
    else stopifnot(all.equal(MM, m, tol = tol))
}

chk.matrix <- function(M) {
    ## check object; including coercion to "matrix" :
    cl <- class(M)
    cat("class ", dQuote(cl), " [",nrow(M)," x ",ncol(M),"]; slots (",
    paste(slotNames(M), collapse=","), ")\n", sep='')
    stopifnot(validObject(M),
          dim(M) == c(nrow(M), ncol(M)),
          identical(dim(m <- as(M, "matrix")), dim(M))
          )
}

isOrthogonal <- function(x, tol = 1e-15) {
    all.equal(diag(as(zapsmall(crossprod(x)), "diagonalMatrix")),
              rep(1, ncol(x)), tol = tol)
}


### ------- Part III --  "Matrix" (classes) specific ----------------------

asD <- function(m) { ## as "Dense"
    if(canCoerce(m, "denseMatrix")) as(m, "denseMatrix")
    else if(canCoerce(m, (cl <- paste(.M.kind(m), "denseMatrix", sep=''))))
        as(m, cl)
    else if(canCoerce(m, "dgeMatrix")) as(m, "dgeMatrix")
    else stop("cannot coerce to a typical dense Matrix")
}

Qidentical <- function(x,y) {
    ## quasi-identical - for 'Matrix' matrices
    if(class(x) != class(y)) return(FALSE)
    slts <- slotNames(x)
    if("factors" %in% slts) { ## allow one empty and one non-empty 'factors'
        slts <- slts[slts != "factors"]
        ## if both are not empty, they must be the same:
        if(length(xf <- x@factors) && length(yf <- y@factors))
            if(!identical(xf, yf)) return(FALSE)
    }
    for(sl in slts)
        if(!identical(slot(x,sl), slot(y,sl)))
            return(FALSE)
    TRUE
}

## Useful Matrix constructors for testing:

rspMat <- function(n, m = n, density = 1/4, nnz = round(density * n*m))
{
    ## Purpose: random sparse Matrix
    ## ----------------------------------------------------------------------
    ## Arguments: (n,m) : dimension [default m=n ==> *square* matrix}
    ##           density: the desired sparseness density:
    ## ----------------------------------------------------------------------
    ## Author: Martin Maechler, Date: 5 Mar 2008, 11:07
    stopifnot(length(n) == 1, n == as.integer(n),
              length(m) == 1, m == as.integer(m),
              0 <= density, density <= 1,
              0 <= nnz, nnz <= n*m)
    x <- numeric(n*m)
    ## entries 2 : (nnz+1) {so we can have '1' as 'special'}
    x[sample(n*m, nnz, replace=FALSE)] <- as.numeric(1L + seq_len(nnz))
    Matrix(x, n,m, sparse=TRUE)
}



rUnitTri <- function(n, upper = TRUE, ...)
{
    ## Purpose: random unit-triangular sparse Matrix .. built from rspMat()
    ## ----------------------------------------------------------------------
    ## Arguments:  n: matrix dimension
    ##         upper: logical indicating if upper or lower triangular
    ##         ...  : further arguments passed to rspMat(), eg. 'density'
    ## ----------------------------------------------------------------------
    ## Author: Martin Maechler, Date:  5 Mar 2008, 11:35

    r <- (if(upper) triu else tril)(rspMat(n, ...))
    ## make sure the diagonal is empty
    diag(r) <- 0
    r <- drop0(r)
    r@diag <- "U"
    r
}

## This is related to rUnitTri(), ver
mkLDL <- function(n, density = 1/3) {
    ## Purpose: make nice artifical   A = L D L'  (with exact numbers) decomp
    ## ----------------------------------------------------------------------
    ## Author: Martin Maechler, Date: 15 Mar 2008
    stopifnot(n == round(n))
    n <- as.integer(n)
    L <- Matrix(0, n,n)

    nnz <- round(n*n * density)
    L[sample(n*n, nnz)] <- seq_len(nnz)
    L <- tril(L,-1)
    diag(L) <- 1
    d.half <- sample(10*(n:1))# random permutation ; use '10*' to be "different" from L entries
    D <- Diagonal(x = d.half * d.half)
    A <- tcrossprod(L * rep(d.half, each=n))
    ## = as(L %*% D %*% t(L), "symmetricMatrix")
    list(A = A, L = L, d.half = d.half, D = D)
}