tests/testthat/test_linalg_closures.R

In mcmcsae: Markov Chain Monte Carlo Small Area Estimation

context("Closures for matrix algebra")

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

test_that("(sparse) matrix sum template works: fixed M0 and variable M1", {
  n <- 10
  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_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)))
})

test_that("(sparse) matrix sum template works: variable M1 and M2", {
  n <- 11
  M1 <- crossprod(rsparsematrix(n, n, density=0.2))
  M2 <- crossprod(rsparsematrix(n, n, density=0.01))
  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_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_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))
  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 <- 11
  M0 <- Diagonal(n)
  M1 <- crossprod(rsparsematrix(n, n, density=0.2))
  M2 <- crossprod(rsparsematrix(n, n, density=0.01))
  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.01))
  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_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_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_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)))
})