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, 0L, 24L, 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=-3L,to=3L, by=0.5)/10
a <- 0.4
b <- 9

tau <- 24L
N.rep <- 1000L
N.boot <- 1000L
outputBase<- "simulation"
                    
Nvals <- c(1000L, 2000L, 5000L, 10000L, 20000L)
Nevents <- as.integer((0.1 + 0.05)*Nvals/2)

dat <- rbind(
             data.frame(varName="", p0=0.1, p1=0.05,
                        beta0=-0.1, beta1=-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)))),
             data.frame(varName="", p0=0.1, p1=0.05,
                        beta0=rep(beta.seq, each=length(beta.seq)),
                        beta1=rep.int(beta.seq, times=length(beta.seq)),
                        N=5000L, Pi=0.04, nBoot=N.boot, nRep=N.rep,
                        event=FALSE),
             data.frame(varName="", p0=0.90, p1=0.95, beta0=-0.1, beta1=-0.1,
                        N=1000L, Pi=0.89, nBoot=N.boot, nRep=N.rep, event=FALSE)
             )

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

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 <- file.path(getwd(), 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=''))

dat<- dat

ncpus <- nrow(dat)

cat(mapply(FUN=function(dat, i) {
  save(dat, a, b, tau,
       file=file.path(new.dir,
         sprintf("%s_input_params_slave%s.rda", outputBase, i)))
#  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/slave_runner.R --vanilla --no-readline --quiet --no-save --args ", outputBase, " ", i,
               " ", new.dir, " 2>&1 | tee ",
               file.path(new.dir,
                         sprintf("%s_output_slave%s.Rout", outputBase, i)),
               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=168:00:00
# You must specify Wall Clock time (hh:mm:ss) [Maximum allowed 30 days = 720:00:00]
#PBS -o %s
# 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.
', file.path(new.dir,
             sprintf("%s_output_cpu%s.output", outputBase, i)),
               cmd)
  pbsFile <- file.path(new.dir,
                       sprintf("%s_output_slave%s.pbs", outputBase, i))
  cat(pbs, file=pbsFile)
##  system(paste("qsub ", pbsFile, sep=''))
  
  cat(cmd, "\n", sep='', append=TRUE, file=output)
  cmd
},
       dat=split(dat, rep(seq_len(ncpus), length.out=nrow(dat))),
       i = sprintf("%0*d", nchar(ncpus, "width"), seq_len(ncpus))), sep='\n')

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)