tests/testthat/test_loglikR.R

In magi: MAnifold-Constrained Gaussian Process Inference

testthat::context("HMC-functions")
library(testthat)
library(magi)

#### set up: taken from HMC-v3 ####
nobs <- 41
noise <- 0.05

VRtrue <- read.csv(system.file("testdata/FN.csv", package="magi"))
pram.true <- list(
  abc=c(0.2,0.2,3),
  rphi=c(0.9486433, 3.2682434),
  vphi=c(1.9840824, 1.1185157),
  sigma=noise
)
fn.true <- VRtrue
fn.true$time <- seq(0,20,0.05)
fn.sim <- fn.true

set.seed(123)
fn.sim[,1:2] <- fn.sim[,1:2]+rnorm(length(unlist(fn.sim[,1:2])), sd=noise)
fn.sim <- fn.sim[seq(1,nrow(fn.sim), length=nobs),]

tvec.nobs <- fn.sim$time
foo <- outer(tvec.nobs, t(tvec.nobs),'-')[,1,]
r <- abs(foo)
r2 <- r^2
signr <- -sign(foo)

expected_marlikmap_par <- c(1.94957912838359, 1.06073179913222, 0.703951749929781, 2.74496621777115,
                            0.0497280144187114)
testthat::test_that("phisigllikC runs without error and is correct", {
  skip_on_os('solaris')  # optim on solaris could fail to find global optimum
  magi:::phisigllikC( c(1.9840824, 1.1185157, 0.9486433, 3.2682434, noise), data.matrix(fn.sim[,1:2]), r)
  fn <- function(par) -magi:::phisigllikC( par, data.matrix(fn.sim[,1:2]), r)$value
  gr <- function(par) -as.vector(magi:::phisigllikC( par, data.matrix(fn.sim[,1:2]), r)$grad)
  marlikmap <- optim(rep(1,5), fn, gr, method="L-BFGS-B", lower = 0.0001)
  testthat::expect_equal(marlikmap$par,
                         expected_marlikmap_par,
                         tolerance = 1e-5)
  
  fn <- function(par) -magi:::phisigllikC( par, data.matrix(fn.sim[,1:2]), r, "compact1")$value
  gr <- function(par) -as.vector(magi:::phisigllikC( par, data.matrix(fn.sim[,1:2]), r, "compact1")$grad)
  marlikmapCompact1 <- optim(rep(1,5), fn, gr, method="L-BFGS-B", lower = 0.0001)
  testthat::expect_equal(marlikmapCompact1$par,
                         c(2.04871398302633, 3.59648132314111, 0.625313733996474, 8.96656240950113, 
                           0.0431289806093459),
                         tolerance = 1e-3)
  
  fn <- function(par) -magi:::phisigllikC( par, data.matrix(fn.sim[,1:2]), r, "rbf")$value
  gr <- function(par) -as.vector(magi:::phisigllikC( par, data.matrix(fn.sim[,1:2]), r, "rbf")$grad)
  marlikmapRbf <- optim(rep(1,5), fn, gr, method="L-BFGS-B", lower = 0.0001)
  testthat::expect_equal(marlikmapRbf$par,
                         c(1.58950017432859, 0.594486181493354, 0.452955008218075, 1.57490985649202, 
                           0.0554703243636422),
                         tolerance = 1e-5)
})


curCovVcompact1 <- calCov(expected_marlikmap_par[1:2], r, signr, kerneltype = "compact1")
curCovRcompact1 <- calCov(expected_marlikmap_par[3:4], r, signr, kerneltype = "compact1")

curCovVrbf <- calCov(expected_marlikmap_par[1:2], r, signr, kerneltype = "rbf")
curCovRrbf <- calCov(expected_marlikmap_par[3:4], r, signr, kerneltype = "rbf")

curCovV <- calCov(expected_marlikmap_par[1:2], r, signr)
curCovR <- calCov(expected_marlikmap_par[3:4], r, signr)

testthat::test_that("calCov runs without error and is correct", {
  varnames <- c("C", "Cprime", "Cdoubleprime", "Cinv", "mphi", "Kphi", "Kinv")
  curCovV.checksum <- sapply(curCovV[varnames], function(x) sum(abs(x)))
  outExpect <- c(387.258476932602, 289.658301140539, 369.687853445841, 
                           1859.64654809337, 224.085486080059, 35.0849556047908, 565.445767848267)
  expect_equal(curCovV.checksum/outExpect, rep(1, length(outExpect)),
               tolerance = 1e-4, check.attributes = FALSE)
  curCovR.checksum <- sapply(curCovR[varnames], function(x) sum(abs(x)))
  outExpect <- c(335.611085502683, 96.4763661864788, 45.1452675162526, 
                 443546.019796851, 244.910143264427, 0.154503453415642, 168699.947543871)
  expect_equal(curCovR.checksum/outExpect, rep(1, length(outExpect)),
               tolerance = 1e-4, check.attributes = FALSE)
  
  curCovV.checksum <- sapply(curCovVcompact1[varnames], function(x) sum(abs(x)))
  outExpect <- c(115.084361649859, 205.278332612955, 2131.40178197008, 
                 37.6922719587419, 144.786133517796, 2075.23766151204, 2.71838147604252)
  expect_equal(curCovV.checksum/outExpect, rep(1, length(outExpect)),
               tolerance = 1e-4, check.attributes = FALSE)
  curCovR.checksum <- sapply(curCovRcompact1[varnames], function(x) sum(abs(x)))
  outExpect <- c(102.781796009223, 103.966939672376, 187.641691435036, 
                 1951.74861141087, 179.878679502857, 42.9071355730163, 120.568332217125)
  expect_equal(curCovR.checksum/outExpect, rep(1, length(outExpect)),
               tolerance = 1e-4, check.attributes = FALSE)
  
  curCovV.checksum <- sapply(curCovVrbf[varnames], function(x) sum(abs(x)))
  outExpect <- c(407.840065042279, 293.008616795784, 338.429857457948, 
                 272835661.009244, 736.757284051975, 0.00565329191451157, 374288748.912434)
  expect_equal(curCovV.checksum/outExpect, rep(1, length(outExpect)),
               tolerance = 1e-4, check.attributes = FALSE)
  curCovR.checksum <- sapply(curCovRrbf[varnames], function(x) sum(abs(x)))
  outExpect <- c(354.860844103142, 95.7538793879521, 43.1253728981135, 
                 786566415.192592, 152.053562220004, 5.00881002256369e-05, 708528973.947655)
  expect_equal(curCovR.checksum/outExpect, rep(1, length(outExpect)),
               tolerance = 1e-4, check.attributes = FALSE)
})
cursigma <- expected_marlikmap_par[5]

testthat::test_that("getMeanCurve runs without error and is correct", {
  startVR <- rbind(magi:::getMeanCurve(fn.sim$time, fn.sim$Vtrue, fn.sim$time,
                                t(expected_marlikmap_par[1:2]), sigma.mat=matrix(cursigma)),
                    magi:::getMeanCurve(fn.sim$time, fn.sim$Rtrue, fn.sim$time,
                                t(expected_marlikmap_par[3:4]), sigma.mat=matrix(cursigma)))
  testthat::expect_equal(sum(startVR), 16.6934654159305, tolerance = 1e-5)
})

testthat::test_that("getMeanCurve deriv is correct", {
  delta <- 0.1
  timeDense <- seq(0, 20, delta)
  ydy <- magi:::getMeanCurve(fn.sim$time, fn.sim$Vtrue, timeDense,
               t(expected_marlikmap_par[1:2]), sigma.mat=matrix(cursigma),
               kerneltype = "generalMatern", deriv = TRUE)
  y <- ydy[[1]]
  dy <- ydy[[2]]
  dyNum <- (y[-(1:2)] - y[-(length(y) - 0:1)]) / (2*delta)
  dy <- dy[c(-1, -length(dy))]
  testthat::expect_equal(dyNum, dy, tolerance = 0.1)
})


testthat::test_that("loglikOrig and loglik runs without error and is correct", {
  out <- magi:::loglikOrig(fn.true[seq(1,nrow(fn.true), length=nobs),],
                    pram.true$abc,
                    c(pram.true$vphi, pram.true$rphi),
                    pram.true$sigma,
                    fn.sim,
                    r,
                    signr)
  expect_equal(out, 
               structure(454.809052651206, 
                         components = structure(c(104.147795635414, 
                                                  108.143924612949, 18.4183735575236, 
                                                  140.020977657482, -1.5461761617613, 
                                                  85.6241573495981), .Dim = 2:3)),
               tolerance = 1e-5)
})

testthat::context("xthetallikC")

dataInput <- data.matrix(fn.sim[,1:2])
xthInit <- c(data.matrix(fn.true[seq(1,nrow(fn.true), length=nobs),1:2]), pram.true$abc)

#### compact1 kernel ####
outExpectedvalue <- -55.73911

testthat::test_that("compact1 - xthetallikC runs without error and is correct", {
  out <- magi:::xthetallikC(dataInput, curCovVcompact1, curCovRcompact1, cursigma, xthInit)
  
  testthat::expect_equal(out$value, outExpectedvalue, tolerance = 1e-4, scale = abs(outExpectedvalue))
  gradExpect <- 143.910609069213
  testthat::expect_equal(sum(out$grad), gradExpect, tolerance = 1e-4, scale = abs(gradExpect))
})

bandsize <- 15
testthat::test_that("compact1 - examine band matrix approximation", {
  curCovVband <- magi:::bandCov(curCovVcompact1, bandsize)
  curCovRband <- magi:::bandCov(curCovRcompact1, bandsize)
  outBandApprox <- magi:::xthetallikBandApproxC(dataInput, curCovVband, curCovRband, cursigma, xthInit)
  testthat::expect_lt(abs((outBandApprox$value - outExpectedvalue)/outExpectedvalue), 1e-3)
  
  delta <- 1e-8
  gradNum <- c()
  for(it in 1:length(xthInit)){
    xthInit1 <- xthInit
    xthInit1[it] <- xthInit1[it] + delta
    gradNum[it] <- 
      (magi:::xthetallikBandApproxC(dataInput, curCovVband, curCovRband, cursigma, xthInit1)$value -
        magi:::xthetallikBandApproxC(dataInput, curCovVband, curCovRband, cursigma, xthInit)$value)/delta
  }
  x <- (gradNum - outBandApprox$grad)/abs(outBandApprox$grad)
  testthat::expect_true(all(abs(x) < 1e-3))
})

#### rbf kernel #### doesn't fit well
outExpectedvalue <- -6017094

testthat::test_that("rbf - xthetallikC runs without error and is correct", {
  out <- magi:::xthetallikC(dataInput, curCovVrbf, curCovRrbf, cursigma, xthInit)
  
  testthat::expect_equal(out$value, outExpectedvalue, tolerance = 1e-4, scale = abs(outExpectedvalue))
  testthat::expect_equal(sum(out$grad), 2781668, tolerance = 1e-1)
})

bandsize <- 15
testthat::test_that("rbf - examine band matrix approximation", {
  curCovVband <- magi:::bandCov(curCovVrbf, bandsize)
  curCovRband <- magi:::bandCov(curCovRrbf, bandsize)
  outBandApprox <- magi:::xthetallikBandApproxC(dataInput, curCovVband, curCovRband, cursigma, xthInit)
  testthat::expect_gt(abs((outBandApprox$value - outExpectedvalue)/outExpectedvalue), 0.25)
  
  delta <- 1e-8
  gradNum <- c()
  for(it in 1:length(xthInit)){
    xthInit1 <- xthInit
    xthInit1[it] <- xthInit1[it] + delta
    gradNum[it] <- 
      (magi:::xthetallikBandApproxC(dataInput, curCovVband, curCovRband, cursigma, xthInit1)$value -
        magi:::xthetallikBandApproxC(dataInput, curCovVband, curCovRband, cursigma, xthInit)$value)/delta
  }
  x <- (gradNum - outBandApprox$grad)/abs(outBandApprox$grad)
  testthat::expect_true(all(abs(x) < 1e-3)) # gradient is self-consistent
})

#### matern kernel ####
outExpectedvalue <- -94.8205825207303

testthat::test_that("xthetallikC runs without error and is correct", {
  out <- magi:::xthetallikC(dataInput, curCovV, curCovR, cursigma, xthInit)
  
  testthat::expect_equal(out$value, outExpectedvalue, tolerance = 1e-5, scale = abs(outExpectedvalue))
  gradExpect <- 167.746373733369
  testthat::expect_equal(sum(out$grad), gradExpect, tolerance = 1e-4, scale = abs(gradExpect))
})

testthat::test_that("xthetallik_rescaledC runs without error and compare to non-scaled", {
  out <- magi:::xthetallik_rescaledC(dataInput, curCovV, curCovR, cursigma, xthInit)
  
  testthat::expect_equal(out$value, outExpectedvalue, tolerance = 1e-5, scale = abs(outExpectedvalue))
  gradExpect <- 167.746373733369
  testthat::expect_equal(sum(out$grad), gradExpect, tolerance = 1e-4, scale = abs(gradExpect))
})

testthat::test_that("xthetallik_rescaledC compare to prior tempered xthetallik", {
  dataInputWithMissing <- dataInput
  dataInputWithMissing[-seq(1,nrow(dataInputWithMissing),4),] <- NA
  out <- magi:::xthetallik_rescaledC(dataInputWithMissing, curCovV, curCovR, cursigma, xthInit)
  pTemp <- nrow(dataInput)/nrow(na.omit(dataInputWithMissing))
  out2 <- magi:::xthetallikC(dataInputWithMissing, curCovV, curCovR, cursigma, xthInit,
                            useBand = FALSE, priorTemperature = c(pTemp, pTemp, 1))
  testthat::expect_equal(out, out2)
})

testthat::test_that("prior tempered xthetallik is the same as rescaling phi", {
  dataInputWithMissing <- dataInput
  dataInputWithMissing[-seq(1,nrow(dataInputWithMissing),4),] <- NA
  
  pTemp <- nrow(dataInput)/nrow(na.omit(dataInputWithMissing))
  out2 <- magi:::xthetallikC(dataInputWithMissing, curCovV, curCovR, cursigma, xthInit,
                            useBand = FALSE, priorTemperature = c(pTemp, pTemp, 1))
  
  
  phiV <- expected_marlikmap_par[1:2]
  phiV[1] <- phiV[1]*pTemp
  phiR <- expected_marlikmap_par[3:4]
  phiR[1] <- phiR[1]*pTemp
  curCovVtempered <- calCov(phiV, r, signr)
  curCovRtempered <- calCov(phiR, r, signr)
  out3 <- magi:::xthetallikC(dataInputWithMissing, curCovVtempered, curCovRtempered, cursigma, xthInit,
                            useBand = FALSE, priorTemperatureInput = 1)
  testthat::expect_equal(out3$value, out2$value, tolerance=1e-3*abs(out2$value))
  testthat::expect_equal(out3$grad, out2$grad, tolerance=0.001)
})

testthat::test_that("prior tempered xthetallik with two separate temperature, and derivatives", {
  dataInputWithMissing <- dataInput
  dataInputWithMissing[-seq(1,nrow(dataInputWithMissing),4),] <- NA
  
  pTemp <- nrow(dataInput)/nrow(na.omit(dataInputWithMissing))
  out <- magi:::xthetallikC(dataInputWithMissing, curCovV, curCovR, cursigma, xthInit,
                           useBand = FALSE, priorTemperature = c(pTemp, pTemp*2))
  out$value
  
  delta <- 1e-8
  gradNum <- c()
  for(it in 1:length(xthInit)){
    xthInit1 <- xthInit
    xthInit1[it] <- xthInit1[it] + delta
    gradNum[it] <- 
      (magi:::xthetallikC(dataInputWithMissing, curCovV, curCovR, cursigma, xthInit1,
                         useBand = FALSE, priorTemperature = c(pTemp, pTemp*2))$value -
         magi:::xthetallikC(dataInputWithMissing, curCovV, curCovR, cursigma, xthInit,
                           useBand = FALSE, priorTemperature = c(pTemp, pTemp*2))$value)/delta
  }
  x <- (gradNum - out$grad)/abs(out$grad)
  testthat::expect_true(all(abs(x) < 1e-3)) # gradient is self-consistent
})

testthat::test_that("xthetallik_withmuC runs without error and compare to zero-mean", {
  out <- magi:::xthetallik_withmuC(dataInput, curCovV, curCovR, cursigma, xthInit)
  out2 <- magi:::xthetallikWithmuBandC(dataInput, curCovV, curCovR, cursigma, xthInit, FALSE)
  testthat::expect_equal(out, out2)
  
  testthat::expect_equal(out$value, outExpectedvalue, tolerance = 1e-5, scale = abs(outExpectedvalue))
  gradExpect <- 167.746373733369
  testthat::expect_equal(sum(out$grad), gradExpect, tolerance = 1e-4, scale = abs(gradExpect))
  
  dotmu <- fnmodelODE(pram.true$abc, data.matrix(fn.true[seq(1,nrow(fn.true), length=nobs),1:2]))
  
  curCovV_withmu <- curCovV
  curCovR_withmu <- curCovR
  curCovV_withmu$mu <- fn.true[seq(1,nrow(fn.true), length=nobs),1]
  curCovR_withmu$mu <- fn.true[seq(1,nrow(fn.true), length=nobs),2]
  curCovV_withmu$dotmu <- dotmu[,1]
  curCovR_withmu$dotmu <- dotmu[,2]
  
  out <- magi:::xthetallik_withmuC(dataInput, curCovV_withmu, curCovR_withmu, cursigma, xthInit)
  out2 <- magi:::xthetallikWithmuBandC(dataInput, curCovV_withmu, curCovR_withmu, cursigma, xthInit, FALSE)
  testthat::expect_equal(out, out2)
  
  dataWN <- dataInput - data.matrix(fn.true[seq(1,nrow(fn.true), length=nobs),1:2])
  xthWN <- c(rep(0,length(xthInit)-3), 0,0,1)
  outWN <- magi:::xthetallikC(dataWN, curCovV, curCovR, cursigma, xthWN)
  
  testthat::expect_equal(out$value, outWN$value, tolerance = 1e-5)
  testthat::expect_equal(out$grad, outWN$grad, tolerance = 1e-5)
})

testthat::test_that("xthetallik_withmuC derivatives", {
  curCovV_withmu <- curCovV
  curCovR_withmu <- curCovR
  curCovV_withmu$mu[] <- mean(dataInput[,"Vtrue"])
  curCovR_withmu$mu[] <- mean(dataInput[,"Rtrue"])
  
  out <- magi:::xthetallik_withmuC(dataInput, curCovV_withmu, curCovR_withmu, cursigma, xthInit)
  out2 <- magi:::xthetallikWithmuBandC(dataInput, curCovV_withmu, curCovR_withmu, cursigma, xthInit, FALSE)
  testthat::expect_equal(out, out2)
  
  delta <- 1e-7
  gradNum <- c()
  for(it in 1:length(xthInit)){
    xthInit1 <- xthInit
    xthInit1[it] <- xthInit1[it] + delta
    gradNum[it] <- 
      (magi:::xthetallik_withmuC(dataInput, curCovV_withmu, curCovR_withmu, cursigma, xthInit1)$value -
        magi:::xthetallik_withmuC(dataInput, curCovV_withmu, curCovR_withmu, cursigma, xthInit)$value)/delta
  }
  x <- (gradNum - out$grad)/abs(out$grad)
  testthat::expect_true(all(abs(x) < 1e-3)) # gradient is self-consistent
  
  dotmu <- fnmodelODE(pram.true$abc, data.matrix(fn.true[seq(1,nrow(fn.true), length=nobs),1:2]))
  
  curCovV_withmu <- curCovV
  curCovR_withmu <- curCovR
  curCovV_withmu$mu <- fn.true[seq(1,nrow(fn.true), length=nobs),1]
  curCovV_withmu$mu <- (curCovV_withmu$mu - mean(curCovV_withmu$mu))*0.5
  curCovR_withmu$mu <- fn.true[seq(1,nrow(fn.true), length=nobs),2]
  curCovR_withmu$mu <- curCovR_withmu$mu - mean(curCovR_withmu$mu)
  curCovV_withmu$dotmu <- dotmu[,1]*0.5
  curCovR_withmu$dotmu <- dotmu[,2]
  
  out <- magi:::xthetallik_withmuC(dataInput, curCovV_withmu, curCovR_withmu, cursigma, xthInit)
  out2 <- magi:::xthetallikWithmuBandC(dataInput, curCovV_withmu, curCovR_withmu, cursigma, xthInit, FALSE)
  testthat::expect_equal(out, out2)
  
  delta <- 1e-7
  gradNum <- c()
  for(it in 1:length(xthInit)){
    xthInit1 <- xthInit
    xthInit1[it] <- xthInit1[it] + delta
    gradNum[it] <- 
      (magi:::xthetallik_withmuC(dataInput, curCovV_withmu, curCovR_withmu, cursigma, xthInit1)$value -
        magi:::xthetallik_withmuC(dataInput, curCovV_withmu, curCovR_withmu, cursigma, xthInit)$value)/delta
  }
  x <- (gradNum - out$grad)
  testthat::expect_true(all(abs(x) < 1e-3)) # gradient is self-consistent
  
})

curCovV_withmu <- curCovV
curCovR_withmu <- curCovR
curCovV_withmu$mu <- sin(fn.sim$time)
curCovV_withmu$mu <- cos(fn.sim$time)
curCovV_withmu$dotmu <- cos(fn.sim$time)
curCovR_withmu$dotmu <- -sin(fn.sim$time)

testthat::test_that("xthetallik_withmuC derivatives - fake mu", {

  out <- magi:::xthetallik_withmuC(dataInput, curCovV_withmu, curCovR_withmu, cursigma, xthInit)
  out2 <- magi:::xthetallikWithmuBandC(dataInput, curCovV_withmu, curCovR_withmu, cursigma, xthInit, FALSE)
  testthat::expect_equal(out, out2)
  
  delta <- 1e-7
  gradNum <- c()
  for(it in 1:length(xthInit)){
    xthInit1 <- xthInit
    xthInit1[it] <- xthInit1[it] + delta
    gradNum[it] <- 
      (magi:::xthetallik_withmuC(dataInput, curCovV_withmu, curCovR_withmu, cursigma, xthInit1)$value -
        magi:::xthetallik_withmuC(dataInput, curCovV_withmu, curCovR_withmu, cursigma, xthInit)$value)/delta
  }
  x <- (gradNum - out$grad)/out$grad
  testthat::expect_true(all(abs(x) < 1e-3)) # gradient is self-consistent
  
  
  out2 <- magi:::xthetallikWithmuBandC(dataInput, curCovV_withmu, curCovR_withmu, cursigma, xthInit, FALSE,
                                priorTemperatureInput = c(4,7))
  
  delta <- 1e-7
  gradNum <- c()
  for(it in 1:length(xthInit)){
    xthInit1 <- xthInit
    xthInit1[it] <- xthInit1[it] + delta
    gradNum[it] <- 
      (magi:::xthetallikWithmuBandC(dataInput, curCovV_withmu, curCovR_withmu, cursigma, xthInit1, FALSE,
                             priorTemperatureInput = c(4,7))$value -
        magi:::xthetallikWithmuBandC(dataInput, curCovV_withmu, curCovR_withmu, cursigma, xthInit, FALSE,
                               priorTemperatureInput = c(4,7))$value)/delta
  }
  x <- (gradNum - out2$grad)/out2$grad
  testthat::expect_true(all(abs(x) < 1e-3)) # gradient is self-consistent
  
})

bandsize <- 10
testthat::test_that("band matrix approximation for withmean", {
  curCovV_withmu <- magi:::bandCov(curCovV_withmu, bandsize)
  curCovR_withmu <- magi:::bandCov(curCovR_withmu, bandsize)
  
  out <- magi:::xthetallik_withmuC(dataInput, curCovV_withmu, curCovR_withmu, cursigma, xthInit)
  out2 <- magi:::xthetallikWithmuBandC(dataInput, curCovV_withmu, curCovR_withmu, cursigma, xthInit, FALSE)
  
  testthat::expect_equal(out, out2)
  
  outWithmeanBand <- magi:::xthetallikWithmuBandC(dataInput, curCovV_withmu, curCovR_withmu, cursigma, xthInit)
  testthat::expect_lt(abs((outWithmeanBand$value - out$value)/out$value), 1e-3)
  
  delta <- 1e-7
  gradNum <- c()
  for(it in 1:length(xthInit)){
    xthInit1 <- xthInit
    xthInit1[it] <- xthInit1[it] + delta
    gradNum[it] <- 
      (magi:::xthetallikWithmuBandC(dataInput, curCovV_withmu, curCovR_withmu, cursigma, xthInit1)$value -
        magi:::xthetallikWithmuBandC(dataInput, curCovV_withmu, curCovR_withmu, cursigma, xthInit)$value)/delta
  }
  x <- (gradNum - outWithmeanBand$grad)/outWithmeanBand$grad
  testthat::expect_true(all(abs(x) < 1e-3)) # gradient is self-consistent
  
  
  outWithmeanBand <- magi:::xthetallikWithmuBandC(dataInput, curCovV_withmu, curCovR_withmu, cursigma, xthInit,
                                           priorTemperatureInput = c(2,3))
  
  delta <- 1e-7
  gradNum <- c()
  for(it in 1:length(xthInit)){
    xthInit1 <- xthInit
    xthInit1[it] <- xthInit1[it] + delta
    gradNum[it] <- 
      (magi:::xthetallikWithmuBandC(dataInput, curCovV_withmu, curCovR_withmu, cursigma, xthInit1,
                             priorTemperatureInput = c(2,3))$value -
        magi:::xthetallikWithmuBandC(dataInput, curCovV_withmu, curCovR_withmu, cursigma, xthInit,
                               priorTemperatureInput = c(2,3))$value)/delta
  }
  x <- (gradNum - outWithmeanBand$grad)/outWithmeanBand$grad
  testthat::expect_true(all(abs(x) < 1e-3)) # gradient is self-consistent
})

bandsize <- 15
curCovVband <- magi:::bandCov(curCovV, bandsize)
curCovRband <- magi:::bandCov(curCovR, bandsize)

testthat::test_that("examine band matrix approximation", {
  outBandApprox <- magi:::xthetallikBandApproxC(dataInput, curCovVband, curCovRband, cursigma, xthInit)
  testthat::expect_lt(abs((outBandApprox$value - outExpectedvalue)/outExpectedvalue), 1e-3)
  
  delta <- 1e-8
  gradNum <- c()
  for(it in 1:length(xthInit)){
    xthInit1 <- xthInit
    xthInit1[it] <- xthInit1[it] + delta
    gradNum[it] <- 
      (magi:::xthetallikBandApproxC(dataInput, curCovVband, curCovRband, cursigma, xthInit1)$value -
        magi:::xthetallikBandApproxC(dataInput, curCovVband, curCovRband, cursigma, xthInit)$value)/delta
  }
  x <- (gradNum - outBandApprox$grad)/abs(outBandApprox$grad)
  testthat::expect_true(all(abs(x) < 1e-3))
})


testthat::test_that("examine low rank approximation", {
  plot(1/curCovV$Keigen1over)
  magi:::truncEigen(1/curCovV$Keigen1over)
  plot(curCovV$KeigenVec[,1])
  plot(curCovV$KeigenVec[,2])
  plot(curCovV$KeigenVec[,3])
  plot(curCovV$KeigenVec[,4])
  
  approxCovVgood <- magi:::truncCovByEigen(curCovV,
                                    magi:::truncEigen(rev(curCovV$Ceigen1over), 0.95),
                                    magi:::truncEigen(rev(curCovV$Keigen1over), 0.95))
  approxCovRgood <- magi:::truncCovByEigen(curCovR,
                                    magi:::truncEigen(rev(curCovR$Ceigen1over), 0.95),
                                    magi:::truncEigen(rev(curCovR$Keigen1over), 0.95))
  
  outApprox <- magi:::xthetallikC(dataInput,
                           approxCovVgood, 
                           approxCovRgood, 
                           cursigma, xthInit)
  testthat::expect_lt((outApprox$value-outExpectedvalue)/outExpectedvalue, 1e-2)
  
  approxCovV <- magi:::truncCovByEigen(curCovV,
                                magi:::truncEigen(rev(curCovV$Ceigen1over), 0.85),
                                magi:::truncEigen(rev(curCovV$Keigen1over), 0.85))
  
  x <- approxCovV$KeigenVec%*%diag(approxCovV$Keigen1over)%*%t(approxCovV$KeigenVec)
  sum((x - approxCovV$Kinv)^2)/sum(approxCovV$Kinv^2)
  diag(abs(x - approxCovV$Kinv)/abs(approxCovV$Kinv))
  
  approxCovR <- magi:::truncCovByEigen(curCovR,
                                magi:::truncEigen(rev(curCovR$Ceigen1over), 0.85),
                                magi:::truncEigen(rev(curCovR$Keigen1over), 0.85))
  
  outApprox <- magi:::xthetallikC(dataInput, approxCovV, approxCovR, cursigma, xthInit)
  abs((outApprox$value - outExpectedvalue)/outExpectedvalue)
  
  delta <- 1e-8
  gradNum <- c()
  for(it in 1:length(xthInit)){
    xthInit1 <- xthInit
    xthInit1[it] <- xthInit1[it] + delta
    gradNum[it] <- 
      (magi:::xthetallikC(dataInput, approxCovV, approxCovR, cursigma, xthInit1)$value -
        magi:::xthetallikC(dataInput, approxCovV, approxCovR, cursigma, xthInit)$value)/delta
  }
  x <- (gradNum - outApprox$grad)/abs(outApprox$grad)
  testthat::expect_true(all(abs(x) < 1e-3))
})



testthat::test_that("band matrix likelihood wrapped runs correctly", {
  datainput <- scan(system.file("testdata/data_band.txt",package="magi"),
                    sep = "\n", what = character(), quiet=TRUE)
  datainput <- strsplit(datainput, "\t")
  datainput <- lapply(datainput, function(x) as.numeric(na.omit(as.numeric(x))))
  
  covVpart <- curCovVband
  covVpart$mphiBand <- matrix(datainput[[2]], nrow=2*datainput[[8]]+1)
  covVpart$KinvBand <- matrix(datainput[[3]], nrow=2*datainput[[8]]+1)
  covVpart$CinvBand <- matrix(datainput[[4]], nrow=2*datainput[[8]]+1)
  covVpart$bandsize <- datainput[[8]]
  
  covRpart <- curCovRband
  covRpart$mphiBand <- matrix(datainput[[5]], nrow=2*datainput[[8]]+1)
  covRpart$KinvBand <- matrix(datainput[[6]], nrow=2*datainput[[8]]+1)
  covRpart$CinvBand <- matrix(datainput[[7]], nrow=2*datainput[[8]]+1)
  covRpart$bandsize <- datainput[[8]]
  yobs <- matrix(datainput[[11]], ncol=2)
  yobs[yobs==-99999] <- NaN
  foo <- magi:::xthetallikBandApproxC(yobs, covVpart, covRpart, datainput[[10]], datainput[[1]])
  
  outsum <- as.numeric(foo$value)+sum(foo$grad)
  testthat::expect_equal(outsum, -655203.16255410481244, tolerance = 1e-3)
})