--- a
+++ b/tests/testthat/test_outbreaker.R
@@ -0,0 +1,442 @@
+context("Test outbreaker")
+current_R_version <- gsub("R version ([0-9.]+) .+$", "\\1", R.version.string)
+
+## test output format ##
+test_that("outbreaker's output have expected format", {
+    ## skip on CRAN
+    skip_on_cran()
+
+
+    ## get data
+    data(fake_outbreak)
+    x <- fake_outbreak
+
+    ## run outbreaker
+    data <- list(dna = x$dna, dates = x$onset, w_dens = x$w)
+    config <- list(n_iter = 10, sample_every = 1, find_import = FALSE)
+    out <- outbreaker(data, config)
+
+
+    out_df <- as.data.frame(out)
+
+    data <- list(dates = x$onset, w_dens = x$w)
+    out2 <- outbreaker(data, config)
+
+    ## check output
+    expect_is(out, "outbreaker_chains")
+    expect_is(out_df, "data.frame")
+    expect_equal(nrow(out), 10)
+    expect_true(!any(is.na(out_df$post)))
+    expect_true(all(out_df$post> -1e30))
+
+})
+
+
+
+
+## test convergence in various settings ##
+test_that("results ok: DNA + time", {
+    ## skip on CRAN
+    skip_on_cran()
+
+
+    ## get data
+    data(fake_outbreak)
+    x <- fake_outbreak
+
+    ## outbreaker DNA + time ##
+    ## analysis
+    set.seed(1)
+    data <- list(dna = x$dna, dates = x$onset, w_dens = x$w)
+    config <- list(n_iter = 5000, sample_every = 50,
+                   init_tree = "seqTrack", find_import = TRUE,
+                   move_pi = FALSE)
+
+    out <- outbreaker(data = data, config = config)
+
+    ## checks
+    out_smry <- summary(out, burnin = 1000)
+
+    ## approx log post values
+    expect_true(min(out_smry$post) > -1250)
+
+    ## at least 80% ancestries correct
+    temp <- mean(out_smry$tree$from==x$ances, na.rm = TRUE)
+    expect_true(temp > .8)
+
+    ## infection date within 3 days on average
+    temp <- mean(abs(out_smry$tree$time - x$onset), na.rm = TRUE)
+    expect_true(temp < 3.5)
+
+    ## mu within reasonable ranges
+    expect_true(out_smry$mu[2] > 1e-5 &&
+                out_smry$mu[5] < 2e-4)
+
+})
+
+
+
+
+
+test_that("results ok: time, no DNA", {
+    ## skip on CRAN
+    skip_on_cran()
+
+
+    ## get data
+    data(fake_outbreak)
+    x <- fake_outbreak
+
+    ## outbreaker time, no DNA ##
+    ## analysis
+    set.seed(1)
+
+    data <- list(dates = x$onset, w_dens = x$w)
+    config <- list(n_iter = 5000, sample_every = 50,
+                   init_tree="star", find_import = FALSE,
+                   move_kappa = FALSE)
+
+    out_no_dna <- outbreaker(data = data, config = config)
+
+    ## checks
+    out_no_dna.smry <- summary(out_no_dna, burnin = 1000)
+
+    ## approx log post values
+    expect_true(min(out_no_dna.smry$post) > -90)
+
+    ## infection date within 3 days on average
+    temp <- mean(abs(out_no_dna.smry$tree$time - x$onset), na.rm = TRUE)
+    expect_true(temp < 3.5)
+
+    ## check that support for ancestries is weak
+    sup <- na.omit(out_no_dna.smry$tree$support)
+    expect_lt(quantile(sup, .9), .55)
+    expect_lt(mean(sup), .35)
+
+})
+
+
+
+test_that("results ok: time, ctd, DNA", {
+    ## skip on CRAN
+    skip_on_cran()
+
+
+    ## get data
+    data(fake_outbreak)
+    x <- fake_outbreak
+
+    ## outbreaker time, ctd, no DNA ##
+    ## analysis
+    ## set the random number generator kind to old R
+    suppressWarnings(RNGversion("3.5.2"))
+    set.seed(1)
+
+    tTree <- data.frame(i = x$ances,
+                        j = 1:length(x$ances))
+
+    ctd <- sim_ctd(tTree, eps = 0.9, lambda = 0.1)
+
+    data <- list(dates = x$onset, w_dens = x$w,
+                 ctd = ctd[,1:2], dna = x$dna)
+    config <- list(n_iter = 5000, sample_every = 50,
+                   init_tree="star", find_import = FALSE,
+                   move_kappa = FALSE)
+
+    out_ctd <- outbreaker(data = data, config = config)
+
+    ## checks
+    out_ctd.smry <- summary(out_ctd, burnin = 1000)
+
+    ## at least 90% ancestries correct
+    temp <- mean(out_ctd.smry$tree$from==x$ances, na.rm = TRUE)
+    expect_true(temp > .9)
+
+    ## eps and lambda parameter estimates within reasonable ranges
+    quant_eps <- quantile(out_ctd$eps, c(0.25, 0.75))
+    expect_true(quant_eps[[1]] > 0.7 &
+                quant_eps[[2]] < 0.95)
+
+    quant_lambda <- quantile(out_ctd$lambda, c(0.25, 0.75))
+    expect_true(quant_lambda[[1]] > 0.05 &
+                quant_lambda[[2]] < 0.35)
+
+    ## approx log post values
+    expect_true(min(out_ctd.smry$post) > -1200)
+
+    ## reset the Random number generator kind
+    RNGversion(current_R_version)
+
+})
+
+
+
+
+test_that("results ok: easy run, no missing cases, no import", {
+    ## skip on CRAN
+    skip_on_cran()
+
+
+    ## get data
+    data(fake_outbreak)
+    x <- fake_outbreak
+
+    ## outbreaker, no missing cases ##
+    ## analysis
+    set.seed(1)
+    config <-  list(n_iter = 5e3, sample_every = 50,
+                    init_tree = "seqTrack", move_kappa = FALSE,
+                    move_pi = FALSE, init_pi = 1,
+                    find_import = FALSE)
+    data <- list(dna = x$dna, dates = x$onset, w_dens = x$w)
+
+    out_no_missing <- outbreaker(data = data, config = config)
+
+    ## checks
+    out_no_missing_smry <- summary(out_no_missing, burnin = 3500)
+
+    ## approx log post values
+    expect_true(min(out_no_missing_smry$post) > -935)
+
+    ## at least 85% ancestries correct
+    temp <- mean(out_no_missing_smry$tree$from==x$ances, na.rm = TRUE)
+    expect_true(temp > .85)
+
+    ## infection datewithin 3 days on average
+    temp <- mean(abs(out_no_missing_smry$tree$time - x$onset), na.rm = TRUE)
+    expect_true(temp < 3.5)
+
+})
+
+
+
+test_that("results ok: no missing cases, detect imported cases",
+{
+    ## skip on CRAN
+    skip_on_cran()
+
+
+    ## get data
+    data(fake_outbreak)
+    x <- fake_outbreak
+
+    ## outbreaker, no missing cases, detect imported cases ##
+    ## analysis
+    set.seed(1)
+    out_with_import <- outbreaker(data = list(dna = x$dna, dates = x$onset,
+                                              w_dens = x$w),
+                                  config = list(n_iter = 5000, sample_every = 100,
+                                                init_tree="star",
+                                                move_kappa = FALSE, move_pi = FALSE,
+                                                init_pi = 1, find_import = TRUE)
+                                  )
+
+    ## checks
+    out_with_import_smry <- summary(out_with_import, burnin = 500)
+    out_with_import_smry$tree$from[is.na(out_with_import_smry$tree$from)] <- 0
+    x$ances[is.na(x$ances)] <- 0
+
+    ## approx log post values
+    expect_true(min(out_with_import_smry$post) > -460)
+
+    ## at least 80% ancestries correct
+    temp <- mean(out_with_import_smry$tree$from==x$ances, na.rm = TRUE)
+    expect_true(temp >= .80)
+
+    ## infection datewithin 3 days on average
+    temp <- mean(abs(out_with_import_smry$tree$time - x$onset), na.rm = TRUE)
+    expect_true(temp < 3.5)
+
+})
+
+
+
+
+
+
+test_that("results ok: kappa and pi", {
+    ## skip on CRAN
+    skip_on_cran()
+
+
+    ## get data
+    onset <- c(0, 2, 6, 13)
+    w <- c(.25, .5, .25)
+
+    ## outbreaker, no missing cases, detect imported cases ##
+    ## analysis
+    set.seed(1)
+
+    data <- list(dates = onset, w_dens = w)
+    config <- list(n_iter = 5000, sample_every = 50, init_tree = "star",
+                   move_kappa = TRUE, move_pi = TRUE, init_pi = 1,
+                   find_import = FALSE, max_kappa = 10)
+
+
+    out <- outbreaker(data, config)
+
+    smry <- summary(out, burnin = 500)
+
+    ## checks
+    expect_equal(smry$tree$from, c(NA, 1, 2, 3))
+    expect_equal(smry$tree$generations, c(NA, 1, 2, 4))
+    expect_true(min(smry$post) > -35)
+    expect_true(all(smry$pi[3:4] > 0.5 & smry$pi[3:4] < 0.8))
+
+})
+
+
+
+
+
+
+test_that("results ok: kappa from genetic data", {
+
+    ## skip on CRAN
+    skip_on_cran()
+
+    dates <- c(0, 5, 15, 20, 30)
+    w <- dnorm(1:10, 5, 1)
+    dna <- matrix(c("A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A",
+                    "T", "T", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A",
+                    "T", "T", "G", "G", "G", "G", "A", "A", "A", "A", "A", "A",
+                    "T", "T", "G", "G", "G", "G", "T", "T", "A", "A", "A", "A",
+                    "T", "T", "G", "G", "G", "G", "T", "T", "G", "G", "G", "G"),
+                  byrow = TRUE, nrow = 5)
+    dna <- ape::as.DNAbin(dna)
+    data <- outbreaker_data(dates = dates, dna = dna, w_dens = w)
+    config <- create_config(init_mu = 2/12, move_mu = FALSE)
+    res <- outbreaker(data, config)
+    expect_equal(summary(res)$tree$generations,
+                 c(NA, 1L, 2L, 1L, 2L))
+
+})
+
+
+
+
+
+
+
+
+
+test_that("results ok: outbreaker with custom priors",
+{
+    ## skip on CRAN
+    skip_on_cran()
+
+
+    ## get data
+    data(fake_outbreak)
+    x <- fake_outbreak
+
+    data <- list(dna = x$dna, dates = x$onset, w_dens = x$w)
+
+    config <-  list(n_iter = 5e3, sample_every = 50,
+                    init_tree = "star", move_kappa = TRUE,
+                    move_pi = TRUE, init_pi = 0.9,
+                    find_import = FALSE)
+
+
+    ## plot(function(x) dbeta(x, 50, 25)) # to see the distribution
+
+    f_pi_1 <- function(x) dbeta(x$pi, 100, 1, log = TRUE) # stronger prior
+    f_pi_2 <- function(x) 0.0 # flat prior
+
+    set.seed(1)
+
+    out_1 <- outbreaker(data, config, priors = list(pi = f_pi_1))
+    out_2 <- outbreaker(data, config, priors = list(pi = f_pi_2))
+
+    smry_1 <- summary(out_1, burnin = 500)
+    smry_2 <- summary(out_2, burnin = 500)
+
+    expect_true(smry_1$pi[2] > 0.9)
+    expect_true(smry_2$pi[2] > 0.6 && smry_2$pi[5] < 0.9)
+
+})
+
+
+
+
+
+
+
+
+test_that("results ok: outbreaker with fixed number returning priors and likelihoods",
+{
+    ## skip on CRAN
+    skip_on_cran()
+
+
+    ## get data
+
+    data(fake_outbreak)
+    x <- fake_outbreak
+
+    data <- list(dna = x$dna, dates = x$onset, w_dens = x$w)
+
+    config <-  list(n_iter = 1000, sample_every = 10,
+                    init_tree = "star", move_kappa = TRUE,
+                    move_pi = TRUE, init_pi = 1,
+                    find_import = FALSE)
+
+
+    ## custom functions
+
+    f1 <- function(x) return(0.0)
+    f2 <- function(x, y) return(0.0)
+    f3 <- function(x) return(1.123)
+
+    priors1 <- custom_priors(mu = f1, pi = f1)
+    priors2 <- custom_priors(mu = f1, pi = f3)
+
+    ll <- custom_likelihoods(genetic = f2,
+                             timing_infections = f2,
+                             timing_sampling = f2,
+                             reporting = f2)
+
+
+    ## tests
+    out1 <- outbreaker(data, config, priors1, ll)
+    out2 <- outbreaker(data, config, priors2, ll)
+
+    expect_true(all(out1$post == 0))
+    expect_true(all(out2$post == 1.123))
+
+})
+
+
+
+
+
+## test consensus tree
+test_that("test consensus trees", {
+
+    ## skip on CRAN
+    skip_on_cran()
+
+
+    ## get data
+    data(fake_outbreak)
+    x <- fake_outbreak
+
+    ## outbreaker DNA + time ##
+    ## analysis
+    set.seed(1)
+    data <- list(dna = x$dna, dates = x$onset, w_dens = x$w)
+    config <- list(n_iter = 5000, sample_every = 50,
+                   init_tree = "seqTrack", find_import = TRUE,
+                   move_pi = FALSE)
+
+    out <- outbreaker(data = data, config = config)
+
+    ## checks
+    out_smry_mpa <- summary(out, burnin = 1000, method = 'mpa')
+    out_smry_dec <- summary(out, burnin = 1000, method = 'decycle')
+
+    ## as there are no cycles, these should be the same
+    testthat::expect_equal(out_smry_mpa$tree$from,
+                           out_smry_dec$tree$from)
+
+})