inst/test-tools.R

In expm: Matrix Exponential, Log, 'etc'

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

source(system.file("test-tools-1.R", package="Matrix"), keep.source=FALSE)

expm.t.identity <- function(x, method,
                            tol = .Machine$double.eps^0.5,
                            check.attributes = FALSE,
                            ...)
{
  ## Purpose: Test the identity   expm(A') = (expm(A))'
  ## ----------------------------------------------------------------------
  ## Arguments: method, ... :          arguments to  expm()
  ##            tol, check.attributes: arguments to  all.equal()
  ## ----------------------------------------------------------------------
  ## Author: Martin Maechler, Date: 23 Feb 2008, 17:26
    ex <- expm::expm(x   , method=method, ...)
    et <- expm::expm(t(x), method=method, ...)
    all.equal(t(ex), et, tolerance = tol, check.attributes = check.attributes)
}


### This is similar to Matrix'  example(spMatrix) :
##' @title random sparse matrix
##' @param nrow,ncol dimension
##' @param ncol
##' @param nnz number of non-zero entries
##' @param density
##' @param rand.x random number generator for 'x' slot
##' @return an  nrow x ncol  matrix
##' @author Martin Maechler, 14.-16. May 2007
rSpMatrix <- function(nrow, ncol = nrow, density, nnz = density*nrow*ncol,
                      sparse = FALSE,
		      rand.x = function(n) round(100 * rnorm(n)))
{
    stopifnot((nnz <- as.integer(nnz)) >= 0,
	      nrow >= 0, ncol >= 0,
	      nnz <= nrow * ncol)
    xx <- rand.x(nnz)
    ## unfortunately, the two resulting matrices might *not* be identical:
    ## because the x's of repeated  (i,j)'s will be *added* for sparse, but not dense:
    ## set.seed(11); m <- rSpMatrix(12, density = 1/10)
    ## set.seed(11); M <- rSpMatrix(12, density = 1/10, sparse=TRUE)
    if(sparse)
	spMatrix(nrow, ncol,
		 i = sample(nrow, nnz, replace = TRUE),
		 j = sample(ncol, nnz, replace = TRUE), x = xx)
    else {
	m <- matrix(0, nrow, ncol)
	m[cbind(i = sample(nrow, nnz, replace = TRUE),
		j = sample(ncol, nnz, replace = TRUE))] <- xx
	m
    }
}

zeroTrace <- function(m) {
    ## Make the {average} trace to 0  -- as it is inside expm(. "Ward77")
    ## This version also works for 'Matrices'
    stopifnot(length(dim(m)) == 2,
              is.numeric(dd <- diag(m)))
    diag(m) <- dd - mean(dd)
    m
}

uniqEntries <- function(m, diagS = FALSE)
{
    ## Purpose: make the non-zero entries of matrix 'm' ``unique''
    ## ----------------------------------------------------------------------
    ## Arguments:
    ## ----------------------------------------------------------------------
    ## Author: Martin Maechler, Date: 26 Feb 2008, 14:40
    m[m > 0] <-  seq_len(sum(m > 0))
    m[m < 0] <- -seq_len(sum(m < 0))
    if(diagS)
        diag(m) <- 10 * sign(diag(m))
    m
}


## This needs "Matrix" package
rMat <- function(n, R_FUN = rnorm,
                 rcondMin = 1.4 * n ^ -1.6226,
                 iterMax = 100)
{
    ## Purpose: random square matrix "not close to singular"
    ## ----------------------------------------------------------------------
    ## Arguments:
    ## NOTE: needs  Matrix::rcond()
    ## ----------------------------------------------------------------------
    ## Author: Martin Maechler, Date: 19 Jan 2008
    ##
    ##--> /u/maechler/R/MM/Pkg-ex/Matrix/rcondition-numb.R  researches rcond( <random mat>)
    ## Result :
    ##   -log[rcond] = log(Kappa) = 1.051 + 1.6226 * log(n)
    ##   ==================================================
    ##   1/rcond = Kappa = exp(1.051 + 1.6226 * log(n))
    ##                   = 2.8605 * n ^ 1.6226
    ##   ==================================================

    ## since we *search* a bit, take a factor ~ 4  higher rcond:
    ##  4 / 2.8605 ~ 1.4 --> default of rcondMin  above

    stopifnot(require("Matrix")) # needs also as(*, ..) etc
    it <- 1
    rcOpt <- 0
    repeat {
        M <- matrix(R_FUN(n^2), n,n)
        if((rc <- Matrix::rcond(M)) >= rcondMin) break
        if(rc > rcOpt) {
            rcOpt <- rc
            M.Opt <- M
        }
        if((it <- it+1) > iterMax) {
            warning("No Matrix found with rcond() >= ",format(rcondMin),
                    "\n Achieved rcond() = ", format(rcOpt),"\n")
            M <- M.Opt
            break
        }
    }
    M
}

doExtras <- interactive() || nzchar(Sys.getenv("R_EXPM_CHECK_EXTRA")) ||
    identical("true", unname(Sys.getenv("R_PKG_CHECKING_doExtras")))