tointegrate<-function(x, beta, tau, a, b) {
  (a/b)*(x/b)^(a-1)*exp(-beta*pmin(x,tau)-(x/b)^a)
}

calc.alpha <- function(p, beta, Pi, tau, a, b) {
  int.val <- integrate(tointegrate, 0, 24, beta=beta, tau=tau, a=a, b=b,
                       subdivisions=2^25, rel.tol=.Machine$double.eps^0.5)$value
  ans <- -log((p-Pi)/(Pi*(int.val + exp(tau * -beta) * (1 - pweibull(tau, a, b)))))
  return(ans)
}

beta.seq <- seq(from=-3,to=3, by=1)/10
a <- 0.4
b <- 9
inf.time<-36

tau <- 24L
N.rep <- 1000L
N.boot <- 1000L
outputBase<- "simulation_followup"
                    
## dat <- rbind(
##              data.frame(p0=0.1, p1=0.05, beta0=c(-0.1), beta1=c(-0.1),
##                         N=c(1000, 2000, 5000, 10000, 20000), Pi=0.04),
##              data.frame(p0=0.1, p1=0.05,
##                         beta0=rep(beta.seq, each=length(beta.seq)),
##                         beta1=rep(beta.seq, times=length(beta.seq)), N=5000, Pi=0.04),
##              data.frame(p0=0.90, p1=0.95, beta0=-0.1, beta1=-0.1, N=1000, Pi=0.89)
##              )

#Nvals <- c(1220, 2440, 6100, 12200)
#Nvals <- c(610, 3050)
Nvals <- c(610, 1220, 2440, 3050, 6100)
#Nevents <- as.integer((0.1 + 0.05)*Nvals/2)
Nevents <- integer(0)
dat <- data.frame(varName=I(""), p0=0.1, p1=0.05, beta0=c(-0.1), beta1=c(-0.1),
                  N=c(Nvals, Nevents), Pi=0.04, nBoot=N.boot, nRep=N.rep,
                  event=rep(c(FALSE, TRUE), times=c(length(Nvals),
                                              length(Nevents))),
                  inf.time=inf.time)
## dat <- data.frame(varName=I(""), p0=0.1, p1=0.05, beta0=c(-0.1), beta1=c(-0.1),
##                   N=c(Nvals, Nevents), Pi=0.04, nBoot=N.boot, nRep=N.rep,
##                   event=rep(c(FALSE, TRUE), times=c(length(Nvals),
##                                               length(Nevents))),
##                   inf.time=inf.time)


rename <- function(x, j) structure(x[,j], names=names(j))

udatc <- unique(rbind(rename(dat, c(p="p0", beta="beta0", Pi="Pi")),
                      rename(dat, c(p="p1", beta="beta1", Pi="Pi"))))

  
alphas <- unlist(do.call(mapply,
                         args=c(list(FUN=function(p, beta, Pi) calc.alpha(p=p, beta=beta, Pi=Pi, tau=tau, a=a, b=b)),
                           udatc, list(SIMPLIFY=FALSE, USE.NAMES=FALSE))))

alphaKey <- sprintf("%a:%a:%a", udatc$p,udatc$beta,udatc$Pi)
dat$alpha0 <- alphas[match(sprintf("%a:%a:%a", dat$p0,dat$beta0,dat$Pi), alphaKey)]
dat$alpha1 <- alphas[match(sprintf("%a:%a:%a", dat$p1,dat$beta1,dat$Pi), alphaKey)]

dat$psi <- with(dat, log(Pi*(1-p0-p1+Pi)/((p0-Pi)*(p1-Pi))))
dat$Pi <- NULL

dat$varName <- with(dat,
                    gsub('-', 'n',
                         paste("results", paste(ifelse(event, "E", "L"), N, sep=''),
                               "b0", beta0, "b1", beta1, "p0", p0, "p1", p1,
                               "Nr", nRep, "Nb", nBoot, sep='_')))

new.dir <- paste("sim", format(Sys.time(), format="%Y-%m-%d-%H-%M-%S"), sep='_')

dir.create(new.dir)

output <- file.path(new.dir, paste(outputBase,"_output_master.txt", sep=''))

ncpus <- nrow(dat)

datChunks <- split(dat, rep(seq_len(ncpus), length.out=nrow(dat)))

mapply(FUN=function(dat, i) {
  save(dat, a, b, tau, file=file.path(new.dir, paste(outputBase, "_input_params_", i, ".rda",sep='')))
#  cmd <- paste("echo \"/usr/bin/R -f slave_runner.R --vanilla --no-readline --quiet --no-save --args ", outputBase, " ", i,
#               " ", new.dir, " 2>&1 | tee ", file.path(new.dir, paste(outputBase, "_output_cpu",i,".Rout", sep='')),"\" | at now", sep='')

  cmd <- paste("R -f /home/dupontct/Bryan/followup_slave_runner.R --vanilla --no-readline --quiet --no-save --args ",
               outputBase, " ", i, " ", new.dir, " 2>&1 | tee ", file.path(new.dir, paste(outputBase, "_output_cpu",i,".Rout", sep='')), sep='')

pbs <- sprintf('#!/bin/sh
# Beginning of PBS batch script.
#PBS -M charles.dupont@vanderbilt.edu
# Status/Progress EMails sent to "my.address@vanderbilt.edu"
#PBS -m bae
# Email generated at b)eginning, a)bort, and e)nd of jobs
#PBS -l nodes=1
# Nodes required (#nodes:#processors per node:CPU type)
#PBS -l mem=512mb
# Total job memory required (specify how many megabytes)
#PBS -l pmem=512mb
# Memory required per processor (specify how many megabytes)
#PBS -l walltime=720:00:00
# You must specify Wall Clock time (hh:mm:ss) [Maximum allowed 30 days = 720:00:00]
#PBS -o cpu%s.output
# Send job stdout to file "myjob.output"
#PBS -j oe
# Send (join) both stderr and stdout to "myjob.output"
cd /home/dupontct/Bryan
%s
# Replace the above echo command with your executable program
# End of PBS batch script.
', i, cmd)
  pbsFile <- file.path(new.dir, paste(outputBase, "_output_cpu", i, ".pbs", sep=''))
  cat(pbs, file=pbsFile)
  system(paste("qsub ", pbsFile, sep=''))
  
  cat(cmd, "\n", sep='', append=TRUE, file=output)
  cmd
},
       dat=datChunks,
       i = seq_along(datChunks))

cat(new.dir, "\n", sep='')
cat("begining calculation run at ", format(Sys.time()), "\n", sep='', file=output, append=TRUE)
cat("\n", append=TRUE, file=output)