library(BatchExperiments)

get_data_wrapper <- function(job, static,..., tr, ip, gr, pon=0, pfi=0){
  stopifnot(all(tr >0, ip>0, gr>tr))
  # It would seem simpler to use the model ~1 insead of ~0 + scale,
  # but only ~0 + scale gives the model matrices rownames (with R <
  # 2.15), which some code depends on.
  get_data(static=static, jump=job$id, trans_rate_model="~0+scale",
           trans_rate_coefs=tr, ip_coefs=ip, growth_rate_coefs=gr,...)
}

reg <- makeExperimentRegistry(id="simStudy",
                              packages=c("EpiGLM", "plyr"))
repls <- 100

nn <- Sys.info()["nodename"]
if (nn == "phylocluster.ccbb.utexas.edu"){
  resources <- list(queue="wilke")
  chunk.size <- ceiling(18 * repls / 24)
} else {
  resources <- list()
  chunk.size <- 1
}

propFreeImp <- list(0, .1, .5)
tmpf <- function(x) {
  r <- data.frame(OutbreakNumber=1, CaseType="person", scale=1,
                  initialInf=1, propOnsetNA=0, propFreeImp=x)
  r$outbreak_code <- paste(r$OutbreakNumber, r$CaseType, sep="_")
  r
}
statics <- lapply(propFreeImp, tmpf)

tmpf <- function(x){
  id <- paste("prob", x$propFreeImp, sep="")
  id <- gsub("\\.", "p", id)
  addProblem(reg, id=id, static=x, dynamic=get_data_wrapper)
  makeDesign(id, exhaustive=list(floor_times=c(TRUE, FALSE),
                   nreps=c(1, 3, 10), tr=0.01, gr=2.1, ip=3,
                   ip_dispersion=.75))
}
prob.designs <- lapply(statics, tmpf)

adf <- data.frame(terms="scale", intercept="no")
addAlgorithm(reg, id="est", fun=get_est)
alg.design <- makeDesign("est", design=adf)

tmpf <- function(x) {
  addExperiments(reg, prob.designs=x, algo.designs=alg.design, repls=repls)
}
lapply(prob.designs, tmpf)

summarizeExperiments(reg)

chunks <- chunk(findExperiments(reg), chunk.size=chunk.size)
submitJobs(reg, ids=chunks, resources=resources)
showStatus(reg)
save.image(file="simStudy.RData")
q('no')

library(EpiGLM)
library(plyr)

tmpf <- function(x) {
  r <- data.frame(OutbreakNumber=1, CaseType=c("pink", "person"), scale=1,
                  initialInf=1, propOnsetNA=0, propFreeImp=x)
  r$outbreak_code <- paste(r$OutbreakNumber, r$CaseType, sep="_")
  r
}

get_data_wrapper <- function(job, static,..., tr, ip, gr, pon=0, pfi=0){
  stopifnot(all(tr >0, ip>0, gr>tr))
  # It would seem simpler to use the model ~1 insead of ~0 + scale,
  # but only ~0 + scale gives the model matrices rownames (with R <
  # 2.15), which some code depends on.
  get_data(static=static, jump=job$id, trans_rate_model="~0+scale",
           trans_rate_coefs=tr, ip_coefs=ip, growth_rate_coefs=gr,...)
}


static <- tmpf(0)
static$OutbreakNumber <- as.integer(static$OutbreakNumber)
dynamic <- get_data_wrapper(job=list(id=1), static=static, nreps=2, tr=.01, ip=1, gr=1)
