tests/testthat/test_mat_sum.R

In mcmcsae: Markov Chain Monte Carlo Small Area Estimation

context("Matrix addition templates")

set.seed(1, kind="Mersenne-Twister", normal.kind="Inversion")

test_that("(sparse) matrix sum template works: fixed M0 and variable M1", {
  n <- 20
  M0 <- Diagonal(n)
  M1 <- Diagonal(x=runif(n))
  ms <- make_mat_sum(M0=M0, M1=M1)
  expect_equivalent(as.matrix(M0 + 4*M1), as.matrix(ms(M1, w1=4)))
  expect_equivalent(as.matrix(M0 + 4*M1), as.matrix(ms(2*M1, w1=2)))
  n <- 12
  M0 <- Diagonal(n)
  M1 <- crossprod(rsparsematrix(n, n, density=0.2))
  ms <- make_mat_sum(M0=M0, M1=M1)
  expect_equivalent(as.matrix(M0 + 4*M1), as.matrix(ms(M1, w1=4)))
  expect_equivalent(as.matrix(M0 + 4*M1), as.matrix(ms(2*M1, w1=2)))
  M0 <- Diagonal(x=rnorm(n))
  ms <- make_mat_sum(M0=M0, M1=M1)
  expect_equivalent(as.matrix(M0 + 4*M1), as.matrix(ms(M1, w1=4)))
  expect_equivalent(as.matrix(M0 + 4*M1), as.matrix(ms(2*M1, w1=2)))
  M0 <- tcrossprod(rep.int(1,n))
  ms <- make_mat_sum(M0=M0, M1=M1)
  expect_is(ms(M1, w1=4), "matrix")
  expect_equivalent(as.matrix(M0 + 4*M1), as.matrix(ms(M1, w1=4)))
  expect_equivalent(as.matrix(M0 + 4*M1), as.matrix(ms(2*M1, w1=2)))
  ms <- make_mat_sum(M0=M0, M1=M1, force.sparse=TRUE)
  expect_is(ms(M1, w1=4), "dsCMatrix")
  expect_equivalent(as.matrix(M0 + 4*M1), as.matrix(ms(M1, w1=4)))
  expect_equivalent(as.matrix(M0 + 4*M1), as.matrix(ms(2*M1, w1=2)))
  M0 <- crossprod(rsparsematrix(n, n, density=0.01))
  ms <- make_mat_sum(M0=M0, M1=M1)
  expect_equivalent(as.matrix(M0 + 4*M1), as.matrix(ms(M1, w1=4)))
  expect_equivalent(as.matrix(M0 + 4*M1), as.matrix(ms(2*M1, w1=2)))
  M1 <- Diagonal(x=rnorm(n))
  M1[3, 3] <- 0
  ms <- make_mat_sum(M0=M0, M1=M1)
  expect_equivalent(as.matrix(M0 + 4*M1), as.matrix(ms(M1, w1=4)))
  expect_equivalent(as.matrix(M0 + 4*M1), as.matrix(ms(2*M1, w1=2)))
})

test_that("(sparse) matrix sum template works: variable M1 and M2", {
  n <- 15
  M1 <- crossprod(rsparsematrix(n, n, density=0.2))
  M2 <- crossprod(rsparsematrix(n, n, density=0.02))
  ms <- make_mat_sum(M1=M1, M2=M2)
  expect_equivalent(as.matrix(2*M1 + 3*M2), as.matrix(ms(M1, M2, w1=2, w2=3)))
  expect_equivalent(as.matrix(2*M1 + 3*M2), as.matrix(ms(4*M1, 3*M2, w1=0.5)))
  M1 <- as.matrix(M1)
  ms <- make_mat_sum(M1=M1, M2=M2)
  expect_is(ms(M1, M2), "matrix")
  expect_equivalent(as.matrix(2*M1 + 3*M2), as.matrix(ms(M1, M2, w1=2, w2=3)))
  expect_equivalent(as.matrix(2*M1 + 3*M2), as.matrix(ms(4*M1, 3*M2, w1=0.5)))
  ms <- make_mat_sum(M1=M1, M2=M2, force.sparse=TRUE)
  expect_is(ms(M1, M2), "dsCMatrix")
  expect_equivalent(as.matrix(2*M1 + 3*M2), as.matrix(ms(M1, M2, w1=2, w2=3)))
  expect_equivalent(as.matrix(2*M1 + 3*M2), as.matrix(ms(4*M1, 3*M2, w1=0.5)))
  M1 <- crossprod(rsparsematrix(n, n, density=0.2))
  M2 <- crossprod(matrix(rnorm(n^2), n))
  ms <- make_mat_sum(M1=M1, M2=M2)
  expect_is(ms(M1, M2), "matrix")
  expect_equivalent(as.matrix(2*M1 + 3*M2), as.matrix(ms(M1, M2, w1=2, w2=3)))
  expect_equivalent(as.matrix(2*M1 + 3*M2), as.matrix(ms(4*M1, 3*M2, w1=0.5)))
  ms <- make_mat_sum(M1=M1, M2=M2, force.sparse=TRUE)
  expect_is(ms(M1, M2), "dsCMatrix")
  expect_equivalent(as.matrix(2*M1 + 3*M2), as.matrix(ms(M1, M2, w1=2, w2=3)))
  expect_equivalent(as.matrix(2*M1 + 3*M2), as.matrix(ms(4*M1, 3*M2, w1=0.5)))
  M2[1, 1] <- M2[5, 4] <- M2[4, 5] <-0
  ms <- make_mat_sum(M1=M1, M2=M2)
  expect_is(ms(M1, M2), "matrix")
  expect_equivalent(as.matrix(2*M1 + 3*M2), as.matrix(ms(M1, M2, w1=2, w2=3)))
  expect_equivalent(as.matrix(2*M1 + 3*M2), as.matrix(ms(4*M1, 3*M2, w1=0.5)))
  ms <- make_mat_sum(M1=M1, M2=M2, force.sparse=TRUE)
  expect_is(ms(M1, M2), "dsCMatrix")
  expect_equivalent(as.matrix(2*M1 + 3*M2), as.matrix(ms(M1, M2, w1=2, w2=3)))
  expect_equivalent(as.matrix(2*M1 + 3*M2), as.matrix(ms(4*M1, 3*M2, w1=0.5)))
  M2 <- Diagonal(x=diag(M2))
  ms <- make_mat_sum(M1=M1, M2=M2)
  expect_equivalent(as.matrix(2*M1 + 3*M2), as.matrix(ms(M1, M2, w1=2, w2=3)))
  expect_equivalent(as.matrix(2*M1 + 3*M2), as.matrix(ms(4*M1, 3*M2, w1=0.5)))
})

test_that("(sparse) matrix sum template works: fixed M0, variable M1 and M2", {
  n <- 20
  M0 <- Diagonal(n)
  M1 <- Diagonal(x=runif(n))
  M2 <- Diagonal(n)
  ms <- make_mat_sum(M0=M0, M1=M1, M2=M2)
  expect_equivalent(as.matrix(M0 + 4*M1 + M2), as.matrix(ms(M1, M2, w1=4)))
  M2 <- as(Diagonal(n), "symmetricMatrix")
  M2[1, 1] <- 0
  ms <- make_mat_sum(M0=M0, M1=M1, M2=M2)
  expect_equivalent(as.matrix(M0 + 4*M1 + M2), as.matrix(ms(M1, M2, w1=4)))
  n <- 18
  M0 <- Diagonal(n)
  M1 <- crossprod(rsparsematrix(n, n, density=0.2))
  M2 <- crossprod(rsparsematrix(n, n, density=0.02))
  ms <- make_mat_sum(M0=M0, M1=M1, M2=M2)
  expect_equivalent(as.matrix(M0 + 2*M1 + 3*M2), as.matrix(ms(M1, M2, w1=2, w2=3)))
  expect_equivalent(as.matrix(M0 + 2*M1 + 3*M2), as.matrix(ms(4*M1, 3*M2, w1=0.5)))
  M0 <- Diagonal(x=rnorm(n))
  ms <- make_mat_sum(M0=M0, M1=M1, M2=M2)
  expect_equivalent(as.matrix(M0 + 2*M1 + 3*M2), as.matrix(ms(M1, M2, w1=2, w2=3)))
  expect_equivalent(as.matrix(M0 + 2*M1 + 3*M2), as.matrix(ms(4*M1, 3*M2, w1=0.5)))
  M0 <- crossprod(rsparsematrix(n, n, density=0.02))
  ms <- make_mat_sum(M0=M0, M1=M1, M2=M2)
  expect_equivalent(as.matrix(M0 + 2*M1 + 3*M2), as.matrix(ms(M1, M2, w1=2, w2=3)))
  expect_equivalent(as.matrix(M0 + 2*M1 + 3*M2), as.matrix(ms(4*M1, 3*M2, w1=0.5)))
  M0 <- tcrossprod(rep.int(1,n))
  ms <- make_mat_sum(M0=M0, M1=M1, M2=M2)
  expect_is(ms(M1, M2), "matrix")
  expect_equivalent(as.matrix(M0 + 2*M1 + 3*M2), as.matrix(ms(M1, M2, w1=2, w2=3)))
  expect_equivalent(as.matrix(M0 + 2*M1 + 3*M2), as.matrix(ms(4*M1, 3*M2, w1=0.5)))
  ms <- make_mat_sum(M0=M0, M1=M1, M2=M2, force.sparse=TRUE)
  expect_is(ms(M1, M2), "dsCMatrix")
  expect_equivalent(as.matrix(M0 + 2*M1 + 3*M2), as.matrix(ms(M1, M2, w1=2, w2=3)))
  expect_equivalent(as.matrix(M0 + 2*M1 + 3*M2), as.matrix(ms(4*M1, 3*M2, w1=0.5)))
  M2 <- Diagonal(x=rnorm(n))
  ms <- make_mat_sum(M0=M0, M1=M1, M2=M2)
  expect_is(ms(M1, M2), "matrix")
  expect_equivalent(as.matrix(M0 + 2*M1 + 3*M2), as.matrix(ms(M1, M2, w1=2, w2=3)))
  expect_equivalent(as.matrix(M0 + 2*M1 + 3*M2), as.matrix(ms(4*M1, 3*M2, w1=0.5)))
  ms <- make_mat_sum(M0=M0, M1=M1, M2=M2, force.sparse=TRUE)
  expect_is(ms(M1, M2), "dsCMatrix")
  expect_equivalent(as.matrix(M0 + 2*M1 + 3*M2), as.matrix(ms(M1, M2, w1=2, w2=3)))
  expect_equivalent(as.matrix(M0 + 2*M1 + 3*M2), as.matrix(ms(4*M1, 3*M2, w1=0.5)))
  ms <- make_mat_sum(M0=M0, M1=M2, M2=M1)
  expect_is(ms(M2, M1, w1=0), "matrix")
  expect_equivalent(as.matrix(M0 + M1), as.matrix(ms(M2, M1, w1=0)))
  expect_equivalent(as.matrix(M0 + 2*M2 + 3*M1), as.matrix(ms(M2, M1, w1=2, w2=3)))
  expect_equivalent(as.matrix(M0 + 2*M2 + 3*M1), as.matrix(ms(4*M2, 3*M1, w1=0.5)))
  ms <- make_mat_sum(M0=M0, M1=M2, M2=M1, force.sparse=TRUE)
  expect_is(ms(M2, M1, w1=0), "dsCMatrix")
  expect_equivalent(as.matrix(M0 + M1), as.matrix(ms(M2, M1, w1=0)))
  expect_equivalent(as.matrix(M0 + 2*M2 + 3*M1), as.matrix(ms(M2, M1, w1=2, w2=3)))
  expect_equivalent(as.matrix(M0 + 2*M2 + 3*M1), as.matrix(ms(4*M2, 3*M1, w1=0.5)))
  M0 <- crossprod(matrix(rnorm(n*n), n, n))
  M0[2, 2] <- 0; M0[4, 4] <- 0
  M1 <- Diagonal(n)
  M1[1, 1] <- 0; M1[3, 3] <- 0
  M2 <- crossprod(rsparsematrix(n, n, 0.1))
  ms <- make_mat_sum(M0=M0, M1=M1, M2=M2)
  expect_is(ms(M1, M2, w2=0), "matrix")
  expect_equivalent(as.matrix(M0 + 4*M1 + M2), as.matrix(ms(M1, M2, w1=4)))
  expect_equivalent(as.matrix(M0), as.matrix(ms(M1, M2, 0, 0)))
  ms <- make_mat_sum(M0=M0, M1=M1, M2=M2, force.sparse=TRUE)
  expect_is(ms(M1, M2, w2=0), "dsCMatrix")
  expect_equivalent(as.matrix(M0 + 4*M1 + M2), as.matrix(ms(M1, M2, w1=4)))
  expect_equivalent(as.matrix(M0), as.matrix(ms(M1, M2, 0, 0)))
})