gpe: Gaussian Process Everything

# fit gp models and check the gp class works as it should

context('model fitting and gp class functions')

# some fake data
n <- 100
df <- data.frame(x = rnorm(n),
                 z = sample(letters[1:10], n, replace = TRUE))
df$y_gaus <- sin(df$x) + rnorm(n, 0, 0.1)
df$y_bern <- rbinom(n, 1, plogis(df$y_gaus))
df$y_pois <- rpois(n, exp(df$y_gaus) * 3)

test_that('basic model fitting works', {
  
  # simple models with all currently available families
  m_gaus <- gp(y_gaus ~ rbf('x') + iid('z') + iid(),
               data = df)
  m_bern <- gp(y_bern ~ rbf('x') + iid('z'),
               data = df,
               family = binomial)
  m_pois <- gp(y_pois ~ rbf('x') + iid('z'),
               data = df,
               family = poisson)
  
})

test_that('model fitting options work', {
  
  # verbosity
  suppressWarnings(m_gaus <- gp(y_gaus ~ rbf('x') + iid('z') + iid(),
                                data = df,
                                verbose = TRUE))
  
  # offsets
  m_gaus <- gp(y_gaus ~ rbf('x') + iid('z') + iid() + offset(x),
               data = df)
  
  # ~~~~~~~~~~~
  # weights
  m_gaus <- gp(y_pois ~ rbf('x') + iid('z') + iid(),
               data = df,
               family = poisson,
               weights = runif(n))

  # wring length
  expect_error(m_gaus <- gp(y_pois ~ rbf('x') + iid('z') + iid(),
                            data = df,
                            family = poisson,
                            weights = runif(3)))
  
  # non-numeric
  expect_error(m_gaus <- gp(y_pois ~ rbf('x') + iid('z') + iid(),
                            data = df,
                            family = poisson,
                            weights = sample(letters, n, replace = TRUE)))

  # non-vector
  expect_error(m_gaus <- gp(y_pois ~ rbf('x') + iid('z') + iid(),
                            data = df,
                            family = poisson,
                            weights = as.list(runif(n))))

  # non finite
  expect_error(m_gaus <- gp(y_pois ~ rbf('x') + iid('z') + iid(),
                            data = df,
                            family = poisson,
                            weights = rep(NA, n)))

  # negative
  expect_error(m_gaus <- gp(y_pois ~ rbf('x') + iid('z') + iid(),
                            data = df,
                            family = poisson,
                            weights = runif(n, -1, 0)))
  
  # ~~~~~~~~~~~
  # mean functions
  

  # not a function
  expect_error(m_gaus <- gp(y_pois ~ rbf('x') + iid('z') + iid(),
                            data = df,
                            family = poisson,
                            mean_function = list()))
  
  # bad function
  expect_error(m_gaus <- gp(y_pois ~ rbf('x') + iid('z') + iid(),
                            data = df,
                            family = poisson,
                            mean_function = function (data) {2 + -3 * data$tomatoes}))

  # bad function
  expect_error(m_gaus <- gp(y_pois ~ rbf('x') + iid('z') + iid(),
                            data = df,
                            family = poisson,
                            mean_function = function (data) {2 + -3 * data$x[1:3]}))
  
  # ~~~~~~
  # check sparse models
  
  # no inducing points provided
  expect_error(m_gaus <- gp(y_gaus ~ rbf('x') + iid('z') + iid(),
                            data = df,
                            inference = 'FITC'))
  
  # now with inducing points
  m_gaus <- gp(y_gaus ~ rbf('x') + iid('z') + iid(),
               data = df,
               inducing_data = df[1:50, ],
               inference = 'FITC')
  
  # ~~~~~
  # check latent inference methods
  
  # currently no Gaussian+Laplace
  expect_error(m_gaus <- gp(y_gaus ~ rbf('x') + iid('z') + iid(),
                            data = df,
                            inducing_data = df[1:50, ],
                            inference = 'Laplace'))
  
  # can't do direct inference on a non-Gaussian likelihood
  expect_error(m_pois <- gp(y_pois ~ rbf('x') + iid('z') + iid(),
                            data = df,
                            family = poisson,
                            inference = 'full'))
  
  # ~~~~~
  # check hyperinference procedures
  
  # these failing due to iid kernel hyperparameter not updatng properly?
  # m_gaus <- gp(y_gaus ~ rbf('x') + iid('z') + iid(),
  #              data = df,
  #              hyperinference = 'BFGS')
  # 
  # m_gaus <- gp(y_gaus ~ rbf('x') + iid('z') + iid(),
  #              data = df,
  #              hyperinference = 'BFGSrestarts')
  
  m_gaus <- gp(y_gaus ~ rbf('x') + iid(),
               data = df,
               hyperinference = 'BFGS')
  
  m_gaus <- gp(y_gaus ~ rbf('x') + iid(),
               data = df,
               hyperinference = 'BFGSrestarts')
  
  m_gaus <- gp(y_gaus ~ rbf('x') + iid('z') + iid(),
               data = df,
               hyperinference = 'none')
  
})

test_that('model prediction works', {
  
  # create some models
  m_gaus <- gp(y_gaus ~ rbf('x') + iid('z') + iid(),
               data = df)
  m_pois <- gp(y_pois ~ rbf('x') + iid('z'),
               data = df,
               family = poisson)
  
  # ~~~~~~~~~
  
  # predict to new df
  pred_gaus_new <- predict(m_gaus, df[1:50, ])
  pred_pois_new <- predict(m_pois, df[1:50, ])
  
  # self predict
  pred_gaus <- predict(m_gaus)
  pred_pois <- predict(m_pois)
  pred_pois_resp <- predict(m_pois, type = 'response')
  pred_pois_resp_sd <- predict(m_pois, type = 'response', sd = TRUE)
  
  # check classes
  expect_equal(class(pred_gaus), 'matrix')
  expect_equal(class(pred_pois), 'matrix')
  expect_equal(class(pred_gaus_new), 'matrix')
  expect_equal(class(pred_pois_new), 'matrix')
  
  # response predictions apparently aren't matrices
  # expect_equal(class(pred_pois_resp), 'matrix')
  
  expect_equal(class(pred_pois_resp_sd), 'list')
  
  # check dimensions
  expect_equal(dim(pred_gaus), c(100, 1))
  expect_equal(dim(pred_pois), c(100, 1))
  expect_equal(dim(pred_gaus_new), c(50, 1))
  expect_equal(dim(pred_pois_new), c(50, 1))
  
})

test_that('model printing works', {
  
  # create some models
  m_gaus <- gp(y_gaus ~ rbf('x') + iid('z') + iid(),
               data = df)
  m_pois <- gp(y_pois ~ rbf('x') + iid('z'),
               data = df,
               family = poisson)

  expected <- c("A Gaussian process model fitted against 100 observations",
                "\tcall: \t\t\t gp(formula = y_gaus ~ rbf(\"x\") + iid(\"z\") + iid(), data = df) ",
                "\tkernel: \t\t(rbf(x) + iid(z)) + iid() ")  
  expect_equal(capture.output(print(m_gaus)),
               expected)
  
})