source("testfunctions.r")
source("test_timecourse_data.r")
wd <- getwd()
setwd("../..")
source("timecourse.r")
source("util.r")
setwd(wd)

test.timecourse.1 <- function() { with.warnings.as.errors(test.timecourse.1.int()) }
test.timecourse.2 <- function() { with.warnings.as.errors(test.timecourse.2.int()) }

test.timecourse.1.int <- function() {
           
  data <- readexpr.func("../../data/endo.txt")
  tme <- c( 0, 0, 0, 0, 2, 2, 2, 2, 4, 4, 4, 6, 6, 6, 9, 9, 9, 9, 24, 24, 24, 24, 0, 0, 0, 0, 2, 2, 2, 2, 4, 4, 4, 4, 6, 6, 6, 6, 9, 9, 9, 9, 24, 24, 24, 24 )
  grp <- c( 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2 )
  ind <- c( 5, 6, 7, 8, 5, 6, 7, 8, 5, 7, 8, 5, 7, 8, 5, 6, 7, 8, 5, 6, 7, 8, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4 )

  stat <<- NULL
  p <<- NULL

  results.timex <<- timex(dat=data$exprs, tme=tme, grp=grp, ind=ind, B=100, null="curve", intercept=FALSE, outfile=NULL, basis="poly", eps=0.05, seed=1, match=NULL)
print("timex finished")  
  p <<- results.timex$p
  stat <<- results.timex$lr
  
  assert(sum(p != expected.p) == 0, "p values not equal to expected")
  assert(sum(abs(stat - expected.stat) > 1e-10) == 0, "stat values not equal to expected")

  results <<- calc.qvalues(p, "timecourse", data$exprs, "smoother", seq(0, 0.95, by=0.05))

  assert(sum(names(results) != c("pi0", "qobj", "q")) == 0, "results list missing some items")

  assert(round(results$pi0,4) == 0.1351)
  assert(sum(abs(results$q - expected.q) > 1e-10) == 0)

  results <<- order.results(stat, p, results$q, data$exprs, data$geneNames, "q", 0.1)

  fname <- "tmpresults_timecourse1.txt"
  expected.fname <- "data/expected_timecourse1.txt"
  savediffexresults.func(fname, results$psig.ord, results$genesig.ord, results$qsig.ord, results$foldc.ord)

  compare.diffex.files(fname, expected.fname)
}

test.timecourse.2.int <- function() {
           
  data <- readexpr.func("data/endo2xfrm.txt")
  tme <- c( 0, 0, 0, 0, 2, 2, 2, 2, 4, 4, 4, 6, 6, 6, 9, 9, 9, 9, 24, 24, 24, 24, 0, 0, 0, 0, 2, 2, 2, 2, 4, 4, 4, 4, 6, 6, 6, 6, 9, 9, 9, 9, 24, 24, 24, 24 )
  grp <- c( 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2 )
  ind <- c( 5, 6, 7, 8, 5, 6, 7, 8, 5, 7, 8, 5, 7, 8, 5, 6, 7, 8, 5, 6, 7, 8, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4 )

  stat <<- NULL
  p <<- NULL

  results.timex <<- timex(dat=data$exprs, tme=tme, grp=grp, ind=ind, B=3, null="curve", intercept=TRUE, outfile=NULL, basis="ncs", eps=0.05, seed=5, match=NULL)
  p <<- results.timex$p
  stat <<- results.timex$lr
  
  results.qvalues <<- calc.qvalues(p, "timecourse", data$exprs, "smoother", seq(0, 0.95, by=0.05))

  assert(sum(names(results.qvalues) != c("pi0", "qobj", "q")) == 0, "results.qvalues list missing some items")

  results.order <<- order.results(stat, p, results.qvalues$q, data$exprs, data$geneNames, "q", 0.1)

  fname <- "tmpresults_timecourse2.txt"
  expected.fname <- "data/expected_timecourse2.txt"
  savediffexresults.func(fname, results.order$psig.ord, results.order$genesig.ord, results.order$qsig.ord, results.order$foldc.ord)

  compare.diffex.files(fname, expected.fname)

  assert(round(results.qvalues$pi0,4) == 0.7859)

  ## Try using a p-value cutoff of 0.00001
  results.order.2 <<- order.results(stat, p, results.qvalues$q, data$exprs, data$geneNames, "q", 0.000001)

  assert(results.order.2$nsig == 0)
}
