tests/testthat/test_phisigllikC.R

In magi: MAnifold-Constrained Gaussian Process Inference

library(magi)
VRtrue <- read.csv(system.file("testdata/FN.csv", package="magi"))

phitrue <- list(
  compact1 = c(2.618, 6.381, 0.152, 9.636),
  rbf = c(0.838, 0.307, 0.202, 0.653),
  matern = c(2.04, 1.313, 0.793, 3.101),
  periodicMatern = c(2.04, 1.313, 9, 0.793, 3.101, 9),
  generalMatern = c(2.04, 1.313, 0.793, 3.101)
)

nobs <- 41
noise <- 0.5

fn.true <- VRtrue
fn.true$time <- seq(0,20,0.05)
fn.sim <- fn.true
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)
signr <- -sign(foo)

for(kerneltype in c("compact1","rbf","matern","periodicMatern","generalMatern")){
  testpoint <- abs(rnorm(length(phitrue[[kerneltype]])+1))
  xc <- magi:::phisigllikC( testpoint, data.matrix(fn.sim[,1:2]), r, kerneltype)
  xr <- magi:::phisigllik( testpoint, fn.sim[,1:2], r, signr, TRUE, kerneltype)
  
  test_that(paste("magi:::phisigllik value -", kerneltype), {
    if(kerneltype=="compact1"){
      skip("compact1 is no longer in use")
    }
    expect_equal(xc$value, as.numeric(xr), tolerance = 1e-4)
    #expect_equal(rep(1, length(xc$grad)), attr(xr, "grad")/as.numeric(xc$grad), tolerance = 0.01)
  })
  
  test_that(paste("magi:::phisigllik gradient -", kerneltype), {
    if(kerneltype=="periodicMatern"){
      skip("periodicMatern magi:::phisigllik gradient has problem for eta, needs debug")
    }
    x0 <- c(phitrue[[kerneltype]]*exp(rnorm(length(phitrue[[kerneltype]]))/10), noise)
    gradTrue <- magi:::phisigllikC(x0, data.matrix(fn.sim[,1:2]), r, kerneltype)$grad
    gradNum <- c()
    for(i in 1:(length(phitrue[[kerneltype]])+1)){
      x1 = x0
      x1[i] = x1[i] + 1e-9
      gradNum[i] <- ((magi:::phisigllikC(x1, data.matrix(fn.sim[,1:2]), r, kerneltype)$value -
                        magi:::phisigllikC(x0, data.matrix(fn.sim[,1:2]), r, kerneltype)$value)/1e-9)
    }
    expect_equal(gradNum, as.numeric(gradTrue), tolerance = 1e-3)
  })
  
  test_that(paste("phisigSample runs without error -", kerneltype), {
    if(kerneltype=="periodicMatern"){
      skip("periodicMatern phisigSample runs with error, needs debug")
    }
    out <- magi:::phisigSample(data.matrix(fn.sim[,1:2]), r, c(phitrue[[kerneltype]], noise),
                 rep(0.03,5), 20, F, kerneltype)
    testthat::expect_equal(length(out$final), 5)
  })
}

for(kerneltype in c("compact1","rbf","matern")){
  testpoint <- abs(rnorm(length(phitrue[[kerneltype]])+1))
  xc <- magi:::phisigllikC( testpoint, data.matrix(fn.sim[,1:2]), r, kerneltype)
  xc2 <- magi:::phisigllikHard2DC( testpoint, data.matrix(fn.sim[,1:2]), r, kerneltype)
  
  test_that(paste("phisigllik returns -", kerneltype), {
    expect_equal(xc, xc2, tolerance = 1e-4)
  })
  
}