library(testthat)
library(RZooRoH)

#test_check("RZooRoH")
context("Running")

test_that("Running one BBB with 1R model ...",{
  my.mod1R <- zoomodel(predefined=FALSE,K=1,krates=c(10))
  freqfile <- (system.file("exdata","typsfrq.txt",package="RZooRoH"))
  typfile <- (system.file("exdata","typs.txt",package="RZooRoH"))
  frq <- read.table(freqfile,header=FALSE)
  bbb <- zoodata(typfile,supcol=4,chrcol=1,poscol=2,allelefreq=frq$V1)
  bbb_results <- zoorun(my.mod1R, bbb, ids = c(1))
  expect_equal(round(bbb_results@modlik[1],1),-6127.2)
  expect_equal(round(bbb_results@realized[1,1],3),0.059)
  expect_equal(round(bbb_results@krates[1,1],1),73.3)
  expect_equal(round(bbb_results@mixc[1,1],3),0.058)
  expect_equal(dim(bbb_results@hbdseg)[1],12)
})

test_that("Running two BBB with 1R model and Nelder-Mead ...",{
  my.mod1R <- zoomodel(predefined=FALSE,K=1,krates=c(10))
  freqfile <- (system.file("exdata","typsfrq.txt",package="RZooRoH"))
  typfile <- (system.file("exdata","typs.txt",package="RZooRoH"))
  frq <- read.table(freqfile,header=FALSE)
  bbb <- zoodata(typfile,supcol=4,chrcol=1,poscol=2,allelefreq=frq$V1)
  bbb_results2 <- zoorun(my.mod1R, bbb, ids = c(1,2), optim_method = "Nelder-Mead",
                         maxiter=1000)
  expect_equal(round(bbb_results2@modlik[1],1),-6127.2)
  expect_equal(round(bbb_results2@modlik[2],1),-5767.7)
  expect_equal(round(bbb_results2@realized[1,1],3),0.058)
  expect_equal(round(bbb_results2@krates[1,1],1),73.1)
  expect_equal(dim(bbb_results2@hbdseg)[1],17)
})

test_that("Running with local HBD probabilities ...",{
  my.mod1R <- zoomodel(predefined=FALSE,K=1,krates=c(10))
  freqfile <- (system.file("exdata","typsfrq.txt",package="RZooRoH"))
  typfile <- (system.file("exdata","typs.txt",package="RZooRoH"))
  frq <- read.table(freqfile,header=FALSE)
  bbb <- zoodata(typfile,supcol=4,chrcol=1,poscol=2,allelefreq=frq$V1)
  bbb_results3 <- zoorun(my.mod1R, bbb, localhbd = TRUE)
  expect_equal(round(min(bbb_results3@hbdp[[1]][1,1:100]),3),0)
  expect_equal(round(min(bbb_results3@hbdp[[2]][1,1:100]),3),0)
  expect_equal(round(min(bbb_results3@hbdp[[3]][1,1:100]),3),0)
  expect_equal(round(min(bbb_results3@hbdp[[4]][1,1:100]),3),0)
  expect_equal(round(min(bbb_results3@hbdp[[5]][1,1:100]),3),0)
  expect_equal(round(max(bbb_results3@hbdp[[1]][1,1:100]),3),0.81)
  expect_equal(round(max(bbb_results3@hbdp[[2]][1,1:100]),3),0.002)
  expect_equal(round(max(bbb_results3@hbdp[[3]][1,1:100]),3),0.05)
  expect_equal(round(max(bbb_results3@hbdp[[4]][1,1:100]),3),0.057)
  expect_equal(round(max(bbb_results3@hbdp[[5]][1,1:100]),3),0.237)
})

test_that("Running with MixKL model ...",{
  freqfile <- (system.file("exdata","typsfrq.txt",package="RZooRoH"))
  typfile <- (system.file("exdata","typs.txt",package="RZooRoH"))
  frq <- read.table(freqfile,header=FALSE)
  bbb <- zoodata(typfile,supcol=4,chrcol=1,poscol=2,allelefreq=frq$V1)
  Mod4L <- zoomodel(K=4,base_rate=10)
  bbb_mod4l <- zoorun(Mod4L, bbb, localhbd = TRUE)
  expect_equal(round(bbb_mod4l@modlik[1],1),-6127.4)
  expect_equal(round(bbb_mod4l@realized[6,1],3),0.100)
  expect_equal(round(bbb_mod4l@realized[1,2],3),0.052)
  expect_equal(round(bbb_mod4l@realized[3,3],3),0.000)
  expect_equal(round(bbb_mod4l@krates[1,1],1),10)
  expect_equal(round(bbb_mod4l@mixc[1,2],3),0.051)
  expect_equal(round(bbb_mod4l@mixc[2,3],3),0.000)
  expect_equal(dim(bbb_mod4l@hbdseg)[1],48)
  expect_equal(max(bbb_mod4l@hbdseg$number_snp),145)
  expect_equal(max(bbb_mod4l@hbdseg$length),48659487)
  expect_equal(round(min(bbb_mod4l@hbdp[[1]][2,1:100]),3),0)
  expect_equal(round(min(bbb_mod4l@hbdp[[2]][2,1:100]),3),0)
  expect_equal(round(min(bbb_mod4l@hbdp[[3]][2,1:100]),3),0)
  expect_equal(round(min(bbb_mod4l@hbdp[[4]][2,1:100]),3),0)
  expect_equal(round(min(bbb_mod4l@hbdp[[5]][2,1:100]),3),0)
  expect_equal(round(max(bbb_mod4l@hbdp[[1]][2,1:100]),3),0.804)
  expect_equal(round(max(bbb_mod4l@hbdp[[2]][2,1:100]),3),0.009)
  expect_equal(round(max(bbb_mod4l@hbdp[[3]][2,1:100]),3),0.147)
  expect_equal(round(max(bbb_mod4l@hbdp[[4]][2,1:100]),3),0.097)
  expect_equal(round(max(bbb_mod4l@hbdp[[5]][2,1:100]),3),0.213)
})