#### 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) }