### A collection of scripts for dealing with ECD model structures.


##Instead of creating the daCapo object, this iterates through the database
##and pulls the data in one parameter at a time.  Thus, it should avoid
##memory problems that daCapo has with large models.
"dsApplyPrint" <-
function(FUN,...,chains, window, studentModelList = listStudentModels(),
         linkModelList = listLinkModels(),trimRegExp=NULL,
         maxparams=99) {
  FUN <- match.fun(FUN)
  ## Argument defaults
  if (missing(chains)) {
    chains <- getChains()[getCycles()[,"chain_last_cycle"]>0]
  }
  if (missing(window)) {
    maxcycle <- min(getCycles()[chains,"chain_last_cycle"])
    burnin <- 1001
    if (maxcycle < 3000) burnin <- 501
    if (maxcycle < 600) burnin <- 1
    window <- c(burnin,maxcycle)
  }
  ## If trimRegExp is non-null, then remove matching expressions from
  ## student model list and link model list
  if (!is.null(trimRegExp)) {
    trimindex <- grep(trimRegExp,studentModelList)
    if (length(trimindex)>0) {
      studentModelList <- studentModelList [-trimindex]
    }
    trimindex <- grep(trimRegExp,linkModelList)
    if (length(trimindex)>0) {
      linkModelList <- linkModelList [-trimindex]
    }
  }
  ###Student Models
  for (model in studentModelList) {
    cat ("\nStudent Model ",model,"\n")
    paramlist <- listParameterNames(model)
    for (pname in paramlist) {
      cat("Values for parameter ",pname,"\n")
      nparam <- nrow(getParameterLabels(model,pname))
      if (nparam > maxparams) {
        cat("Too many parameters (",nparam,"), skipping.\n")
      } else {
        param <- getParameterVector(model,pname,window,chains)
        print(try(FUN(param,...)))
        cat("\n")
      }
    }
    cat("\n")
  }
  ###Link Models
  for (model in linkModelList) {
    cat ("\nLink Model ",model,"\n")
    paramlist <- listParameterNames(model)
    for (pname in paramlist) {
      cat("Values for parameter ",pname,"\n")
      nparam <- nrow(getParameterLabels(model,pname))
      if (nparam > maxparams) {
        cat("Too many parameters (",nparam,"), skipping.\n")
      } else {
        param <- getParameterVector(model,pname,window,chains)
        print(try(FUN(param,...)))
        cat("\n")
      }
    }
    cat("\n")
  }
}

##Instead of creating the daCapo object, this iterates through the database
##and pulls the data in one parameter at a time.  Thus, it should avoid
##memory problems that daCapo has with large models.
"dsApply" <-
function(SUMMARY, FUN,...,chains, window, studentModelList = listStudentModels(),
         linkModelList = listLinkModels(),trimRegExp=NULL,
         maxparams=99, na.rm=TRUE) {
  ## Leave as string, see docs for do.call (R version 2.0.x)
  ## SUMMARY <- match.fun(SUMMARY)
  FUN <- match.fun(FUN)
  ## Argument defaults
  if (missing(chains)) {
    chains <- getChains()[getCycles()[,"chain_last_cycle"]>0]
  }
  if (missing(window)) {
    maxcycle <- min(getCycles()[chains,"chain_last_cycle"])
    burnin <- 1001
    if (maxcycle < 3000) burnin <- 501
    if (maxcycle < 600) burnin <- 1
    window <- c(burnin,maxcycle)
  }
  ## If trimRegExp is non-null, then remove matching expressions from
  ## student model list and link model list
  if (!is.null(trimRegExp)) {
    trimindex <- grep(trimRegExp,studentModelList)
    if (length(trimindex)>0) {
      studentModelList <- studentModelList [-trimindex]
    }
    trimindex <- grep(trimRegExp,linkModelList)
    if (length(trimindex)>0) {
      linkModelList <- linkModelList [-trimindex]
    }
  }
  modelList <- c(studentModelList,linkModelList)
  ##Loop through models
  results <- lapply(modelList,dsApplyModel,SUMMARY,FUN,...,
                    window=window,chains=chains,maxparams=maxparams,na.rm=na.rm)
  if (na.rm) {
    results <- results[!is.na(results)]
  }
  do.call(SUMMARY,results)
}

## Auxiliary function for doing one model at a time.
"dsApplyModel" <-
function (model, SUMMARY, FUN, ..., window,chains,maxparams=99,na.rm=TRUE) {
  results <- NULL
  paramlist <- listParameterNames(model)
  for (pname in paramlist) {
    nparam <- nrow(getParameterLabels(model,pname))
    if (nparam < maxparams) {
      param <- getParameterVector(model,pname,window,chains)
      if (length(list(...)) > 0) {
        val <- try(FUN(param,...))
      } else {
        ## No optional arguments, call without ... in case function
        ## doesn't accept
        val <- try(FUN(param))
      }
      if (class(val) == "try-error") val <- NA #R
      if (class(val) == "try") val <- NA       # S-Plus
      if (!na.rm || !is.na(val)) {
        if (is.null(results)) {
          results <- list(val)
        } else {
          results <- c(results,val)
        }
      }
    }
  }
  #cat("Results for model ",model,":  ")
  #print(results)
  do.call(SUMMARY,results)
}

## Applies a function recursively to the parts of the daCapo object
## and assumes each result returns a plot. It adds a sub argument to
## the plot based on the model.
## The arguments whichStudentModels and whichLinkModels can be used
## to print out only a subset of the models.
## If prefix is supplied, each model is plotted to a separate PDF file,
## If ask is true, plots are made with
"dsApplyPlot" <-
function(FUN,...,
         chains, window, studentModelList = listStudentModels(),
         linkModelList = listLinkModels(),trimRegExp=NULL,
         maxparams=99, prefix="",ask=missing(prefix),silent=TRUE)
{
  FUN <- match.fun(FUN)
  onscreen <- missing(prefix)
  ## Argument defaults
  if (missing(chains)) {
    chains <- getChains()[getCycles()[,"chain_last_cycle"]>0]
  }
  if (missing(window)) {
    maxcycle <- min(getCycles()[chains,"chain_last_cycle"])
    burnin <- 1001
    if (maxcycle < 3000) burnin <- 501
    if (maxcycle < 600) burnin <- 1
    window <- c(burnin,maxcycle)
  }
  ## If trimRegExp is non-null, then remove matching expressions from
  ## student model list and link model list
  if (!is.null(trimRegExp)) {
    trimindex <- grep(trimRegExp,studentModelList)
    if (length(trimindex)>0) {
      studentModelList <- studentModelList [-trimindex]
    }
    trimindex <- grep(trimRegExp,linkModelList)
    if (length(trimindex)>0) {
      linkModelList <- linkModelList [-trimindex]
    }
  }
  first <- TRUE
  for (model in studentModelList) {
    if (!onscreen) {
      pdf(paste(prefix,model,".pdf",sep=""))
    }
    if (!silent) {
      cat ("\nStudent Model ",model,"\n")
    }
    if (!first && ask) {
      readline("Hit return for next model: ")
    }
    first <- TRUE
    paramlist <- listParameterNames(model)
    for (pname in paramlist) {
      if (!silent) {
        cat("Variable ",pname,"\n")
      }
      nparam = nrow(getParameterLabels(model,pname))
      if (nparam > maxparams) {
        cat("Too many parameters (",nparam,"), skipping variable ",pname," model ",model,"\n")
      } else {
        param <- getParameterVector(model,pname,window,chains)
        if (!first && ask) {
          readline("Hit return for next variable: ")
        }
        first <- FALSE
        try(FUN(param,...,ask=ask,
                sub=paste("Mod: ",model,"Var: ",getObservableVar(pname),sep=" ")))
      }
    }
    if (!onscreen) {
      dev.off()
    }
  }

  for (model in linkModelList) {
    if (!onscreen) {
      pdf(paste(prefix,model,".pdf",sep=""))
    }
    if (!silent) {
      cat ("\nLink Model ",model,"\n")
    }
    if (!first && ask) {
      readline("Hit return for next model: ")
    }
    first <- TRUE
    paramlist <- listParameterNames(model)
    lmname <- splitlmn(model)
    for (pname in paramlist) {
      if (!silent) {
        cat("Variable ",pname,"\n")
      }
      nparam <- nrow(getParameterLabels(model,pname))
      if (nparam > maxparams) {
        cat("Too many parameters (",nparam,"), skipping variable ",pname," model ",model,"\n")
      } else {
        param <- getParameterVector(model,pname,window,chains)
        if (!first && ask) {
          readline("Hit return for next variable: ")
        }
        first <- FALSE
        try(FUN(param,...,ask=ask,
                sub=paste("Task: ",lmname$task,"Var: ",getObservableVar(pname),sep=" ")))
      }
    }
    if (!onscreen) {
      dev.off()
    }
  }

}

### These are functions related to building an "Result Matrix" directly from the database
"dsBuildRMatrix" <-
function(FUN,...,chains, window, 
         linkModelList = listLinkModels(),trimRegExp=NULL,
         maxparams=99,silent=TRUE) {
  FUN <- match.fun(FUN)
  ## Argument defaults
  if (missing(chains)) {
    chains <- getChains()[getCycles()[,"chain_last_cycle"]>0]
  }
  if (missing(window)) {
    maxcycle <- min(getCycles()[chains,"chain_last_cycle"])
    burnin <- 1001
    if (maxcycle < 3000) burnin <- 501
    if (maxcycle < 600) burnin <- 1
    window <- c(burnin,maxcycle)
  }
  ## If trimRegExp is non-null, then remove matching expressions from
  ## student model list and link model list
  if (!is.null(trimRegExp)) {
    trimindex <- grep(trimRegExp,linkModelList)
    if (length(trimindex)>0) {
      linkModelList <- linkModelList [-trimindex]
    }
  }
  lmpnames <- lapply(linkModelList, listParameterNames)
  ## Calculate dimensions of resulting array
  nrow <- do.call("sum",lapply(lmpnames,length))
  irow <- 1
  task <- character(0)
  EM <- character(0)
  observable <- character(0)
  ## We will create the result array when we have the first row calculated. 
  R <- NULL
  for (model in linkModelList) {
    lmname <- splitlmn(model)
    if (!silent) {
      cat("Task: ",lmname$task,"  EM: ",lmname$EM,"\n")
    }
    pnames <- listParameterNames(model)
    task <- c(task,rep(lmname$task,length(pnames)))
    EM <- c(EM,rep(lmname$EM,length(pnames)))
    observable <- c(observable, getObservableVar(pnames))
    lmdat <- getAllParameters(model,window,chains,maxparams)
    results <- t(sapply(lmdat,FUN,...))
    if (is.null(R)) {
      ## Shape results matrix based on first row
      R <- results
    } else {
      R <- rbind(R,results)
    }
  }
  result <- data.frame(task=I(task),EM=I(EM),observable=I(observable),R)
  result
}