#########################################################
## Copyright (c) 2005
## University of Washington
## Licensed under the terms set forth by University of
## Washington. If you did not sign such a license, you
## are using this software/code illegally and you do not
## have permission to use, modify, or redistribute
## this or any files in this software package.
#########################################################

require(tcltk, quietly = TRUE, keep.source = FALSE)  || stop("TCLTK support is absent.")
source("includeui.r")
source("include.r")
source("thread.r")
edgedir <- getwd()

#########################################################################
#       This function computes p-values for temporal differential       #
#         expression.                                                   #
#       Author(s): John Storey, Jeffrey Leek                            #
#       Language: R                                                     #
#       References:                                                     #
#         Storey, J.D., Leek, J.T., Xiao, W., Dai, J.Y., Davis R.W.,    #
#         A Significance Method for Time Course Microarray Experiments  #
#         Applied to Two Human Studies, University of Washington        #
#         Biotatistics Department Technical report (2004),              #
#         http://faculty.washington.edu/~jstorey/                       #
#       Arguments:                                                      #
#         dat - The m x n matrix of expression values.                  #
#         tme - The time covariate vector.                              #
#         grp - If appropriate, the vector of group labels.             #
#         ind - If appropriate, the vector of individual labels.        #
#         B - The number of null permutations.                          #
#         null - The null hypothesis that is to be tested.              #
#         dfo - The degrees of freedom for the spline.                  #
#         intercept - Specifies if an intercept is to be included.      #
#         outfile - If the results should be written to a file (not     #
#           used in the GUI).                                           #
#         basis - The spline basis, either natural cubic spline or      # 
#           polynomial spline.                                          #
#         eps - likelihood based convergence criterion when EM          #
#           algorithm is necessary                                      #
#         seed - option to set seed for random bootstrap samples        #
#                                                                       #
#########################################################################
timex <- function(dat, tme, grp = NULL, ind = NULL, B = 100, 
                  null = c("curve", "linear", "flat"), dfo = NULL,
                  intercept = TRUE, outfile = NULL,
                  basis = c("ncs", "poly"), eps = .05,
                  seed = NULL, match=NULL,
                  updatefunc=NULL) {

  result <- timecourse.func(dat=dat, tme=tme, grp=grp, ind=ind,
                            B=B, null=null, intercept=intercept,
                            outfile=outfile, basis=basis, eps=eps,
                            seed=seed, match=match,
                            background=FALSE)

  return(result)
}

#########################################################################
# This function calls the "EMFitR" C function.
# Arguments:
#########################################################################
em.fit = function(dat, H, ind, eps) {
  m = nrow(dat)
  n = ncol(dat)
  N = length(unique(ind))

  if(!is.loaded("EMFitR"))
    try.load(paste(edgedir, "/edge", .Platform$dynlib.ext, sep = ""))

  out = .C("EMFitR",
    dat = as.double(dat),
    H = as.double(H),
    ind = as.integer(ind),
    idx = as.integer(order(ind) - 1),
    eps = as.double(eps),
    fit = double(m * n),
    alpha = double(m * n),
    m = as.integer(m),
    n = as.integer(n),
    N = as.integer(N),
    nGrp = as.integer(table(ind)),
    iter = integer(1))

  return(list("fit" = matrix(out$fit, nrow = m), "alpha" = matrix(out$alpha, nrow = m), "iter" = out$iter))
}

ps <- function(tme, dfo, intercept = TRUE) {
  xx <- matrix(rep(1, length(tme)), ncol = 1)
  for(i in 1:dfo)
    xx <- cbind(xx, tme ^ i)
  if(!intercept)
    xx <- xx[, -1]
  
  return(xx)
}

#########################################################################
# This function starts a timecourse score calculation in a background thread.
# Arguments:
#   dat - the data to work with
#   updatefunc - a function to call periodically (approx. every 500 ms)
#     allowing a GUI to check up on the thread's status
#   finalfunc - the function to call when the thread completes
#   printmsg - a function that prints its arguments (to console, gui, 
#     whatever
#########################################################################
timecourse.func <- function(
                            dat, tme, grp = NULL, ind = NULL, B = 100, 
                            null = c("curve", "linear", "flat"), dfo = NULL,
                            intercept = FALSE, outfile = NULL,
                            basis = c("ncs", "poly"), eps = .05,
                            seed = NULL,match=NULL,
                            updatefunc=NULL, finalfunc=NULL,
                            printmsg=function(...){},
                            background=TRUE) 
{
  ## If debugging turned on, write the params of this function to a
  ## debug file
  debug.params()

  ## check data
  if(is.null(dat))
    stop("dat is null")
  assign("tdat", dat, .GlobalEnv)

  null <- match.arg(null)
  basis <- match.arg(basis)
# match <- match.arg(match)
  if(!is.null(seed)) {set.seed(seed)}
  require(stats, quietly = TRUE, keep.source = FALSE)
  require(splines, quietly = TRUE, keep.source = FALSE)
  require(mgcv, quietly = TRUE, keep.source = FALSE)
  
  
  
  ## transform these variables into values 1, 2, ...  
  if(!is.null(grp)) {
    grp <- make.consecutive.int(grp)
  }
  if(!is.null(ind)) {
    ind <- make.consecutive.int(ind)
    N <- max(ind)
    if(N < 2) {
      stop("You must have more than one individual under longitudinal sampling.")
    }
  }
  if(!is.null(match)) {
    match <- make.consecutive.int(match)
  }
  
  ## get dimensions
  m <- nrow(dat)
  n <- ncol(dat) #=length(tme)=length(grp)=length(ind)
  if(is.null(grp)) {
    grp <- rep(1, n)
  }
  g <- max(grp)
  
  if(null!="curve" && g>1 && is.null(match))
    stop("Must set null=curve when there are two or more classes.")  
#  if(null!="curve" && g>1) {
#    printmsg(c("Warning: null=", null, "used...ignoring group assignment."), bell=TRUE)
#  }
  if(null=="curve" && g==1)
    stop("Must set null to flat or linear for a within class analysis.")
  
  # Check on matching
  if(!is.null(match)) {
    if(g!=2) { 
      stop("Cannot do a matched analysis with other than 2 classes.")
    }
    if((n%%2) !=0 ) {
      stop("Number of arrays is not even and you have designated a matched design.")
    }
    if(abs(1-min(sort(match[grp==1])==sort(match[grp==2])))) {
      stop("The matching variable is not correctly coded.")
    }
    if(!is.null(ind)) {
      for(i in 1:max(match)){
        if(abs(1-min(sort(tme[grp==1][match[grp==1]==i])==sort(tme[grp==2][match[grp==2]==i])))) {
          stop("The time points on each individual must coincide for a matched design.")
        }
      }
    }  else    {
      if(abs(1-min(sort(tme[grp==1])==sort(tme[grp==2])))) {
        stop("The time points must be the same across classes for a matched design.")
      }
    }

    index <- 1
    xx <- 0*dat[,1:(n/2)]
    for(i in 1:max(match)) {
      matchsize <- sum(match[grp==2]==i)
      xx[,index:(index+matchsize-1)] <- dat[,grp==2][,match[grp==2]==i] - dat[,grp==1][,match[grp==1]==i]
      index <- index + matchsize
    }
    dat <- xx
    rm(xx)
    n <- n/2
    ind <- ind[grp==1]
    tme <- tme[grp==1]
    grp <- rep(1,n)
    ind <- make.consecutive.int(ind)
    N <- max(ind)
    g <- 1
    null <- "flat"
  }


  ## This checks to make sure that there are at least two time points per group
    
  ## Automatically choose the basis dimension for the spline
  ## Ignores intercept, so that doesn't matter at this point in the function

dfo1 <- dfo
  if(is.null(dfo1)) {
    if(is.null(ind)){
      vv <- diag(rep(1, g))[, grp]
      if(g==1)
        vv <- matrix(vv, nrow = 1)
      dat.bar <- t((1 / apply(vv, 1, sum)) * (vv %*% t(dat)))
      dat1 <- dat - dat.bar[, grp]
      for(j in 1:g) {
	if(length(unique(tme[grp==j])) < 2 ){
          stop("There must be at least two unique time points in each group for a time course analysis.");
        }
        else if(length(unique(tme[grp==j])) == 2){
          dfo1 <- 1; 
          break; 
        }
        else if(length(unique(tme[grp==j])) == 3){
          dfo1 <-2; 
          next;
        }
        else{
          oo <- svd(dat1[, grp == j])
          for(i in 1:min(5,sum(grp==j)))
            dfo1 <- max(floor(smooth.spline(tme[grp==j], oo$v[, i])$df-1), dfo1)
        }
      }
      rm(dat1)
    }
    else {
      for(k in 1:N){
        if(length(unique(tme[ind==k])) < 2){
          stop("There must be at least two unique time points for each individual.");
        }
      }
      vv <- diag(rep(1, N))[, ind]
      dat.bar <- t((1 / apply(vv, 1, sum)) * (vv %*% t(dat)))
      dat1 <- dat - dat.bar[, ind]
      
      for(j in 1:g) {
        if(length(unique(tme[grp==j])) < 2 ){
          stop("There must be at least two unique time points in each group",
               "for a time course analysis.");
        }
        else if(length(unique(tme[grp==j])) == 2){
          dfo1 <- 1; 
          break; 
        }
        else if(length(unique(tme[grp==j])) == 3){
          dfo1 <-2; 
          next;
        }
        else{
          oo <- svd(dat1[, grp == j])
          for(i in 1:min(5,sum(grp==j)))
            dfo1 <- max(floor(smooth.spline(tme[grp==j], oo$v[, i])$df-1), dfo1)
        }
      }
      rm(dat1)
    }
  }


  ## Check to make sure that basis dimension is within the proper range
  for(i in 1:g) {
    if(length(unique(tme[grp == i])) == length(tme[grp == i])) {
       dfo1 <- min(length(unique(tme[grp == i])) - 2, dfo1)
     } else {
       dfo1 <- min(length(unique(tme[grp == i])) - 1, dfo1)
     }
  }
  dfo1 <- max(dfo1,1)
  
  if(!is.null(dfo) && (dfo1 != dfo))
  {
    printmsg(paste("Dimension of basis",
                   dfo,
                   "was outside the allowable bounds. It was reset to",
                   dfo1, "."))
  } else {
    printmsg(paste("Dimension of basis chosen:", dfo1), bell=FALSE)
  }
  dfo <- dfo1

  if(g == 1 && !is.null(ind)) {
    intercept <- FALSE
  }
  if(is.null(ind) && null=="flat") {
    intercept <- TRUE
  }
  
  ## set up basis matrices
  ## alternative model
  if(basis == "ncs") {
    knts <- quantile(tme, probs = seq(0, 1, length = (dfo + 1)))[-c(1, (dfo + 1))]
    xx <- ns(tme, knots = knts, intercept = intercept)
  } else if(basis == "poly") {
    xx <- ps(tme, dfo = dfo, intercept = intercept)
  }
  ## null model
  xx0 <- xx
  if(null == "flat") {
    xx0 <- matrix(rep(1, n), ncol = 1)
  } 
  if(null == "linear") {
    xx0 <- ns(tme, intercept = intercept)
  }

  ### get fitted values, residual sum of squares, F-stats
  ### independent sampling
  if(is.null(ind)) {

    if(!intercept) {    
      vv <- diag(rep(1, g))[, grp]
      if(g == 1) 
        vv <- matrix(vv, nrow = 1)
      dat.bar <- t((1 / apply(vv, 1, sum)) * (vv %*% t(dat)))
      dat <- dat - dat.bar[, grp]
      for(i in 1:g) {
        for(j in 1:ncol(xx)) {
          xx[grp == i, j] <- xx[grp == i, j] - mean(xx[grp == i, j])
        }
      }
      for(i in 1:g) {
        for(j in 1:ncol(xx0)) {
          xx0[grp == i, j] <- xx0[grp == i, j] - mean(xx0[grp == i, j])
        }
      }
    }
    
    if(intercept || (null != "flat")) {
      H0 <- xx0 %*% solve(t(xx0) %*% xx0) %*% t(xx0)
    } else {
      H0 <- matrix(rep(0, n * n), ncol = n)
    }
    
    H1 <- 0 * H0
    for(i in 1:g) {
      H1[grp == i, grp == i] <- xx[grp == i, ] %*% solve(t(xx[grp == i, ]) %*% xx[grp == i, ]) %*% t(xx[grp == i, ])
    }
    
    fit1 <- t(H1 %*% t(dat))
    res = dat - fit1
    rss1 = drop((res ^ 2) %*% rep(1, n))
    fit0 <- t(H0 %*% t(dat))
    rss0 = drop(((dat - fit0) ^ 2) %*% rep(1, n))
  }
  
  ##########################################################################

  ### get fitted values, residual sum of squares, F-stats
  ### longitudinal sampling
  if(!is.null(ind)) {
    if(intercept) {
      H0 <- xx0 %*% solve(t(xx0) %*% xx0) %*% t(xx0)
      H1 <- 0 * H0
      for(i in 1:g)
        H1[grp == i, grp == i] <- xx[grp == i, ] %*% solve(t(xx[grp == i, ]) %*% xx[grp == i, ]) %*% t(xx[grp == i, ])
      ef0 <- em.fit(dat, H0, ind, eps = eps)
      ef1 <- em.fit(dat, H1, ind, eps = eps)
      fit0 <- ef0$fit 
      alpha0 <- ef0$alpha
      fit1 <- ef1$fit 
      alpha1 <- ef1$alpha
      res0 <- dat - fit0 - alpha0
      rss0 = drop(((dat - fit0) ^ 2) %*% rep(1, n))
      rss1 = drop(((dat - fit1) ^ 2) %*% rep(1, n))
    }
    if(!intercept) {
      vv <- diag(rep(1, N))[, ind]
      if(is.vector(vv)) vv <- matrix(vv, byrow=TRUE, nrow=1)
      dat.bar <- t((1 / apply(vv, 1, sum)) * (vv %*% t(dat)))
      dat <- dat - dat.bar[, ind]
      for(i in 1:N) {
        for(j in 1:ncol(xx)) {
          xx[ind == i, j] <- xx[ind == i, j] - mean(xx[ind == i, j])
        }
      }
      H0 <- matrix(rep(0, n * n), ncol = n)
      if(null != "flat") {
        for(i in 1:N) {
          for(j in 1:ncol(xx0)) {
            xx0[ind == i, j] <- xx0[ind == i, j] - mean(xx0[ind == i, j])
          }
        }
        H0 <- xx0 %*% solve(t(xx0) %*% xx0) %*% t(xx0)
      }
      H1 <- 0 * H0
      
      for(i in 1:g) {
        H1[grp == i, grp == i] <- xx[grp == i, ] %*% solve(t(xx[grp == i, ]) %*% xx[grp == i, ]) %*% t(xx[grp == i, ])
      }
    
      fit1 <- t(H1 %*% t(dat))
      res <- dat - fit1
      rss1 = drop((res ^ 2) %*% rep(1, n))
      fit0 <- t(H0 %*% t(dat))
      rss0 = drop(((dat - fit0) ^ 2) %*% rep(1, n))

      w <- diag(length(tme))
      for(i in 1:N)
        w[ind == i, ind == i] <- w[ind == i, ind == i] - 1 / sum(ind == i)
      
      rmv <- rep(0, N)
      for(i in 1:N)
        rmv[i] <- which.max(ind == i)
# Take away length(rmv) columns (ie length(rmv)*nrow(dat) data values)
      res0 <- res[, -rmv]
      wsq <- mat.sq(w[-rmv, -rmv])
      wsq.inv <- solve(wsq)
      res00 <- res0 %*% wsq.inv
    }
  }

  FF <- (rss0 - rss1) / rss1

  ptm <- proc.time()[3]
  lr0 <- NULL
  tcks <- round(c(B/3, 2*B/3))

  Sys.sleep(0.1)
  
  if(!is.null(seed)) {set.seed(seed)}

  if(is.null(ind)) 
  # Create a n * B matrix
    v = matrix(sample(1:n, n * B, replace = TRUE), nrow = B)
  else if(!intercept)
  # Create a (n-N)*B matrix
    v = matrix(sample(1:ncol(res00), ncol(res00) * B, replace = TRUE), nrow = B)
  else {
    v = sample(1:n, 2 * n * B, replace = TRUE)
    v1 = v[((n * B) + 1):(2 * n * B)]
    v = matrix(v[1:(n * B)], nrow = B)
    v1 = matrix(v1, nrow = B)
  }
 
# gah - move m, n up here to get into datalist 
  m <- nrow(dat)
  n <- ncol(dat)

  datalist <- list()
if(debugOn) print("begin datalist creation")
  datalist$B <- as.integer(B)
  datalist["updatefunc"] <- list(updatefunc)
  datalist["finalfunc"] <- list(finalfunc)
  datalist["printmsg"] <- list(printmsg)
  datalist$m <- as.integer(m)
  datalist$n <- as.integer(n)
  datalist$dat <- dat
  datalist$startTime = proc.time()[3]
  datalist["intercept"] <- list(intercept)

# gah - the following notation ensures that null 
# elements get added to the list

  if(exists("ind") && !is.null(ind)) datalist["ind"] <- list(as.integer(ind))   
  else datalist["ind"] <- list(NULL)
  datalist["fit0"] <- list(fit0)
  datalist["v"] <- list(as.integer(v-1))
  datalist["H0"] <- list(H0)
  datalist["H1"] <- list(H1)
  datalist["FF"] <- list(FF)

# gah - the following are either added or null is added to the
# list to be sure that position numbers stay constant.

  if(exists("res") && !is.null(res)) datalist$res <- res
  else datalist["res"] <- list(NULL)
  if(exists("res0") && !is.null(res0)) datalist$res0 <- res0
  else datalist["res0"] <- list(NULL)
  if(exists("res00") && !is.null(res00)) datalist$res00 <-  res00
  else datalist["res00"] <- list(NULL)
  if(exists("rmv") && !is.null(rmv)) datalist$rmv <- as.integer(rmv-1)
  else datalist["rmv"] <- list(NULL)
  if(exists("alpha1") && !is.null(alpha1)) datalist$alpha1 <- alpha1
  else datalist["alpha1"] <- list(NULL)
  if(exists("eps") && !is.null(eps)) datalist$eps <- eps
  else datalist["eps"] <- list(NULL)
  if(exists("wsq") && !is.null(wsq)) datalist$wsq <- wsq
  else datalist["wsq"] <- list(NULL)
  if(exists("N") && !is.null(N)) datalist$N <- as.integer(N)
  else datalist["N"] <- list(NULL)
  if(exists("v1") && !is.null(v1)) datalist$v1 <- as.integer(v1-1) # n*B
  else datalist["v1"] <- list(NULL)
  if(exists("rmv")) datalist["inc"] <- list(as.integer((1:n)[-rmv]-1))
  else datalist["inc"] <- list(NULL)
if(debugOn) print("add Idx")
  if(exists("ind") && !is.null(ind)) {
     datalist["Idx"] <- list(as.integer(order(ind)-1))
     datalist["nGrp"] <- list(as.integer(table(ind)))
  }
  else {
     datalist["Idx"] <- list(NULL)
     datalist["nGrp"] <- list(NULL)
  }

# gah - debug datalist
if(debugOn) {
   print("print datalist")
   if(debugOn==2) print(datalist)
  for(i in 1:length(datalist)) {
     if(is.vector(datalist[[i]])) print(paste(i,
          typeof(datalist[[i]]),names(datalist[i]),
          length(datalist[[i]]),datalist[[i]][[1]]))
     else print(paste(i,typeof(datalist[[i]]),names(datalist[i]),
          length(datalist[[i]])))
      }
   }
if(debugOn) print(length(datalist))
if(debugOn) savedata <<- datalist;

  if (background)
    {
      check.not.nullna(datalist$updatefunc)

      ## Start a background thread using timecourse.start()
      ## Note that endfn=timecourse.return.func means that timecourse.return.func()
      ## will be called when the thread finishes.
      ##
      ## Check the thread's progress every 500 ms (1/2 second),
      ## calling updatefunc at each check.
      tid = thread(timecourse.start, startdata=datalist, updatefn=datalist$updatefunc,
        updatedata=datalist, endfn=timecourse.return.func, enddata=datalist,
        updateInterval=500)

      return(tid)
    }
  else
    {
      result <- timecourse.start(datalist, background=FALSE)

      pvals <- get.pvalues(datalist$FF, result$F)
      return(list(lr=datalist$FF, lr0=result$F, p=pvals))
    }
  
}

#########################################################################
# This function calls the low-level C function RunTimecourse.
# Arguments:
#   datalist - a list with data needed by the C function.
#########################################################################
timecourse.start <- function(datalist, background=TRUE)
{
  if(!is.list(datalist))
    stop("datalist argument [passed to timecourse.start()] is not a list")

  dat <- check.not.nullna(datalist$dat)

  m = nrow(dat)
  n = ncol(dat)

  ind <- datalist$ind
  intercept <- datalist$intercept
  if(!is.logical(intercept))
    stop("could not find datalist$intercept (or it is not a logical value)")

  if(!is.loaded("RunTimecourseIndep"))
    try.load(paste(edgedir, "/edge", .Platform$dynlib.ext, sep = ""))

  if(is.null(ind)) 
  {
    # User chose "None (independent sampling)" under "Choose covariate corresponding to individuals"
#    tid <- .Call("RunTimecourseIndep", 
#      check.not.nullna(datalist$fit0),
#      check.not.nullna(datalist$res),
#      check.not.nullna(datalist$v)-1,
#      check.not.nullna(datalist$H1),
#      check.not.nullna(datalist$H0),
#      check.not.nullna(datalist$B),
#      m, n, background)      
# gah - change to pass datalist instead of separate arguments

# gah - check for null or na arguments, but don't use the return value
# also, this is a reminder of which list elements will be used.
      check.not.nullna(datalist$fit0)
      check.not.nullna(datalist$res)
      check.not.nullna(datalist$v)
      check.not.nullna(datalist$H1)
      check.not.nullna(datalist$H0)
      check.not.nullna(datalist$B)
    tid <- .Call("RunTimecourseIndep", datalist, background)
  }
  else if(!intercept)
  {
    # User did NOT check "Check to include baseline levels in analysis" 
    check.not.nullna(datalist$rmv)
#    tid <- .Call("RunTimecourseNoInt",
#      fit = check.not.nullna(datalist$fit0),
#      res = check.not.nullna(datalist$res00),
#      v = check.not.nullna(datalist$v)-1,
#      rmv = rmv -1,
#      inc = (1:n)[-rmv] - 1,
#      H1 = check.not.nullna(datalist$H1),
#      H0 = check.not.nullna(datalist$H0),
#      wsq = check.not.nullna(datalist$wsq),
#      ind = ind,
#      idx = order(ind) - 1,
#      m = m,
#      n = n,
#      N = check.not.nullna(datalist$N),
#      nGrp = table(ind),
#      check.not.nullna(datalist$B),
#      background)
# gah - change to pass datalist instead of separate arguments
# gah - check for null or na arguments, but don't use the return value
# also, this is a reminder of which list elements will be used.
      check.not.nullna(datalist$fit0)
      check.not.nullna(datalist$res00)
      check.not.nullna(datalist$v)
      check.not.nullna(datalist$H1)
      check.not.nullna(datalist$H0)
      check.not.nullna(datalist$wsq)
      check.not.nullna(datalist$N)
      check.not.nullna(datalist$B)
    tid <- .Call("RunTimecourseNoInt", datalist, background)
  }
  else
  {
    # User DID check "Check to include baseline levels in analysis"  
#    tid <- .Call("RunTimecourseInt",
#      fit = check.not.nullna(datalist$fit0),
#      res = check.not.nullna(datalist$res0),
#      v = check.not.nullna(datalist$v)-1,
#      v1 = check.not.nullna(datalist$v1)-1,
#      H1 = check.not.nullna(datalist$H1),
#      H0 = check.not.nullna(datalist$H0),
#      alpha1 = check.not.nullna(datalist$alpha1),
#      ind = ind,
#      idx = order(ind) - 1,
#      eps = check.not.nullna(datalist$eps),
#      m = m,
#      n = n,
#      N = check.not.nullna(datalist$N),
#      nGrp = table(ind),
#      B = check.not.nullna(datalist$B),
#      background)

# gah - change to pass datalist instead of separate arguments
# gah - check for null or na arguments, but don't use the return value
# also, this is a reminder of which list elements will be used.
      check.not.nullna(datalist$fit0)
      check.not.nullna(datalist$res0)
      check.not.nullna(datalist$v)
      check.not.nullna(datalist$v1)
      check.not.nullna(datalist$H1)
      check.not.nullna(datalist$H0)
      check.not.nullna(datalist$alpha1)
      check.not.nullna(datalist$eps)
      check.not.nullna(datalist$N)
      check.not.nullna(datalist$B)

    tid <- .Call("RunTimecourseInt", datalist, background)

    }

  return(tid)
}

#########################################################################
# This function is used internally. It is called when the background
# timecourse thread finishes.
#########################################################################
timecourse.return.func <- function(tid, progress, datalist)
{
  if(as.logical(progress$canceled))
    return("canceled")

  if(!is.list(datalist))
    stop("datalist passed to timecourse.return.func is not a list")

  # Unwrap datalist.
  finalfunc <- check.not.nullna(datalist$finalfunc)
  printmsg <- check.not.nullna(datalist$printmsg)
  
  FF <- check.not.nullna(datalist$FF)
  result <- check.not.nullna(timecourse.getresult(tid))

  F.0 <- result$F
  if(is.null(F.0))
    stop("result list did not contain a member named F")

  progress$succeeded=TRUE
  
  ## NaN or NA results invalidate the whole analysis
  if(any(is.na(FF)) || any(is.na(F.0))) {
    print("Found NaN/NA...")
    progress$succeeded=FALSE
    progress$reason=paste("The analysis encountered a divide-by-zero error.",
      "This probably means there is a systematic error in your data resulting",
      "in a zero variance for a gene.")
    finalfunc(p=NULL, lr=NULL, lr0=NULL, progress)
    return(0)
  }

  pvals <- get.pvalues(FF, F.0)
  finalfunc(p=pvals, lr=FF, lr0=F.0, progress)
}

timecourse.getresult <- function(tid)
{
  lis <- .Call("GetTimecourseResult", tid=tid)
  return(lis)
}

