tests/matrix-ex.R

In Rmpfr: R MPFR - Multiple Precision Floating-Point Reliable

stopifnot(require("Rmpfr"))

x <- mpfr(0:7, 64)/7
mx <- x
dim(mx) <- c(4,2)
(m. <- mx) # "print"
m.[,2] <- Const("pi", 80)
m.[,] <- exp(mpfr(1, 90))
stopifnot(is(mx, "mpfrMatrix"), dim(mx) == c(4,2),
          is(m., "mpfrMatrix"), dim(m.) == dim(mx),
	  dim(is.finite(mx)) == dim(mx),
	  dim(is.nan(mx)) == dim(mx),
          getPrec(m.) == 90)

xx <- (0:7)/7
m.x <- matrix(xx, 4,2)
m2 <- mpfr(xx, 64); dim(m2) <- dim(m.x)
##
u <- 10*(1:4)
y <- 7 * mpfr(1:12, 80)
my <- y
dim(my) <- 3:4
m.y <- asNumeric(my)
stopifnot(all.equal(m2, mpfr(m.x, 64), tol=0), # not identical(..)
	  my[2,2] == 35,
	  my[,1] == 7*(1:3))

.N <- function(x) { if(!is.null(dim(x))) as(x,"array") else as(x,"numeric") }
noDN <- function(.) { dimnames(.) <- NULL ; . }
allEQ <- function(x,y) all.equal(x,y, tol=1e-15)

## FIXME write "functions" that take  x -> {mx , m.x}  and run the following as *function*
## ----  then use previous case *and* cases with NA's !
##       and use  higher precision via  fPrec = 2 etc ...

stopifnot(allEQ(m.x, noDN(.N(mx))),
	  allEQ(m.y, noDN(.N(my))),
	  allEQ(noDN(.N(my %*% mx)), m.y %*% m.x),
	  allEQ(noDN(.N(crossprod(mx, t(my)))), crossprod(m.x, t(m.y))),
	  allEQ(noDN(.N(tcrossprod(my, t(mx)))),
			tcrossprod(m.y, t(m.x))),
	  ##
	  identical(mx, t(t(mx))),
	  identical(my, t(t(my))),
	  ## matrix o vector .. even  vector o vector
	  identical(noDN(.N(my %*% 1:4)), m.y %*% 1:4 ),
	  identical(noDN(.N(my %*% my[2,])), m.y %*% .N(my[2,])),
	  identical( crossprod(1:3, my), 1:3 %*%   my),
	  identical(tcrossprod(1:4, my), 1:4 %*% t(my)),
	  identical(crossprod(y), t(y) %*% y),
	  identical(tcrossprod(y), y %*% t(y)),
	  identical(noDN(.N( crossprod(y))), crossprod(7 * 1:12)),
	  identical(noDN(.N(tcrossprod(y))),tcrossprod(7 * 1:12)),
	  identical(tcrossprod(1:3, u), noDN(.N(tcrossprod(1:3, as(u,"mpfr")))))
	  )

mx[3,1] <- Const("pi", 64)
stopifnot(allEQ(sum(mx[,1]), pi + 4/7))
m2 <- mx[c(1,4),]
stopifnot(dim(m2) == c(2,2), sum(m2) == 2)

## "mpfrArray" o "mpfr" :
Tmx <- array(TRUE, dim(mx), dimnames=dimnames(mx))
stopifnot(identical(Tmx, mx == (mx - mpfr(0, 10))),
	  identical(Tmx, mx - mpfr(1, 10) * mx == 0))
## subassignment, many kinds
mx[5] <- pi
mx[6] <- Const("pi",100)
stopifnot(validObject(mx), allEQ(mx[5], mx[6]),
	  getPrec(mx) == c(rep(64,5), 100, 64,64))

## %*% with vectors on LHS, ...
y <- t(2:4) # 1 x 3 matrix
m1 <-     (0:10)     %*% y
m2 <- mpfr(0:10, 50) %*% y
stopifnot((d <- m1 - m2) == 0, identical(dim(m1), dim(d)),
          m2 == m1, m1 == m2)

r <- 10*(0:4)
y <- t(2:6)
m1 <- 1:3 %*% y  %*% r
y. <- t(mpfr(2:6, 20))
m2 <- 1:3 %*% y. %*% r
stopifnot(m1 == m2, m1 - m2 == 0, identical(dim(m1), dim(m2)))

### Array (non-matrix) ---- indexing & sub-assignment :
A <- mpfrArray(1:24, prec = 96, dim = 2:4)
a <-     array(1:24,            dim = 2:4)
a.1 <- as(A[,,1], "array")
a1. <- as(A[1,,], "array")
A1. <- as(A[1,,], "mpfr")

stopifnot(all.equal(noDN(a.1), a[,,1], tol=0),
	  identical(A1., as.vector(A[1,,])),
	  ## arithmetic, subsetting etc:
	  allEQ(noDN(.N(A / A1.)), a/c(a1.)),
	  allEQ(noDN(.N(a / A1.)), a/c(a1.)),
          identical(noDN(A == 23), a == 23),
          identical(noDN(10 >= A), 10 >= a),
          identical(noDN(A <=  2), a <=  2),
          identical(noDN(A < 2.5), a < 2.5),
          identical(noDN(A !=  5), a !=  5),
          identical(A != 3, !(3 == A)),
          identical(-1 > A, A == 100),
          identical(noDN(A <= 0), a == pi)
	  )

A[1,2,3] <- Const("pi")
A[1, ,2] <- 1 / A[1,,2]
A

## check that A is "==" a  where
a <- array(1:24, 2:4); a[1,2,3] <- pi; a[1,,2] <- 1/a[1,,2]
stopifnot(allEQ(noDN(.N(A)), a),
          ## check aperm() methods :
          allEQ(noDN(.N(aperm(A))), aperm(a)),
          {p <- c(3,1:2); allEQ(noDN(.N(aperm(A,p))), aperm(a,p))},
          {p <- c(2:1,3); allEQ(noDN(.N(aperm(A,p))), aperm(a,p))})

## cbind() / rbind():
options(warn = 2)## no warnings here - ("exact recycling"):
validObject(m0 <- cbind(pi=pi, i = 1:6))
validObject(m1 <- cbind(a=Const("pi",60),i = 1:6, "1/mp" = 1/mpfr(1:3,70)))
validObject(m2 <- cbind(pi=pi, i = 1:2, 1/mpfr(1:6,70)))
validObject(n2 <- rbind(pi=pi, i = 1:2, 1/mpfr(1:6,70)))
stopifnot(is(m0,"matrix"), is(m1, "mpfrMatrix"), is(m2, "mpfrMatrix"),
          dim(m0) == c(6,2), dim(m1) == c(6,3), dim(m2) == c(6,3))
options(warn = 1)
suppressWarnings(eval(ex <- quote(m3 <- cbind(I=10, 1:3, inv=1/mpfr(2:3,80)))))
validObject(suppressWarnings(     n3 <- rbind(I=10, 1:3, inv=1/mpfr(2:3,80))))
stopifnot(identical(t(n2), m2),
          identical(t(n3), m3), validObject(m3),
          is(tryCatch(eval(ex), warning=function(.).), "warning"),
	  identical(cbind("A", "c"), matrix(c("A", "c"), 1,2)),
	  identical(rbind("A", 2),   matrix(c("A", "2"), 2,1)) )

## head() / tail() :
stopifnot(all.equal(c(21, 12),
		    dim(mm3 <- m3[rep(1:3, each=7), rep(3:1, 4)])),
	  all.equal(dim(t3 <- tail(mm3)), c(6, 12)),
	  all.equal(head(mm3), mm3[1:6,]))

## matrix(<mpfr>) works since 2015-02-28:
x <- mpfr(pi,64)*mpfr(2,64)^(2^(0:19))
(mx <- matrix(x, 4,5))

stopifnot(is(mx, "mpfrMatrix"),
    all.equal(matrix(0:19, 4,5),
              asNumeric(log2(log2(mx) - log2(Const("pi")))),
              tol = 1e-15)) # 64b-lnx: see 8.1e-17


## Ensure that apply() continues to work with 'bigz'/'bigq'
A <- matrix(2^as.bigz(1:12), 3,4)
Aq <- as.bigq(A)
mA <- as(A, "mpfr") # failed {as dim(A) is "double", not "integer"}
(Qm <- A / (A^2 - 1)) # big rational matrix
MQ <- mpfr(Qm, precBits = 64)
stopifnot(exprs = {
    mA == A
    mA == Aq
    is.bigq(Aq)
    mA == mpfr(A, precBits=16)
    mA == asNumeric(A)
    is.bigq(Qm)
    is(MQ, "mpfrMatrix")
    all.equal(Qm, MQ, tol = 1e-18)
    identical(dim(mA), dim(A))
    identical(dim(mA), dim(Qm))
    identical(asNumeric(apply(A,  1, min)),
              apply(asNumeric(A), 1, min))
    identical(asNumeric(apply(A,  1, max)),
              apply(asNumeric(A), 1, max))
    identical(asNumeric(apply(A,  2, max)),
              apply(asNumeric(A), 2, max))
    identical(asNumeric(apply(A,  2, min)),
              apply(asNumeric(A), 2, min))
})
## mA etc, failed up to Rmpfr 0.8-1; the apply() parts failed up to Rmpfr 0.6.0

if(FALSE) ## Bug in gmp :  apply(*, range) does *not* return matrix
stopifnot(
          identical(asNumeric(apply(A,  1, range)),
                    apply(asNumeric(A), 1, range))
)

if(FALSE) ## Bug in gmp : --- no mean method for bigz, just mean.bigq
stopifnot(
          all.equal(asNumeric(apply(A,  1, mean)),
                    apply(asNumeric(A), 1, mean))
          ,
          all.equal(asNumeric(apply(A,  2, mean)),
                    apply(asNumeric(A), 2, mean))
)

cat('Time elapsed: ', proc.time(),'\n') # "stats"
if(!interactive()) warnings()