forwsim <- function(dta.N=2,th.u=1,th.v=1,th.au=3,th.bu=1,th.av=3,th.bv=1,dsg.uth=0,dsg.hyper=FALSE,dsg.pool=FALSE,dsg.logsd=TRUE,dsg.patlimit=-1,dsg.kini=1,dsg.c0=0,dsg.c1=1,dsg.c2=10000,dsg.p0=0.5,dsg.a=.05,dsg.b=.1,dsg.s0=.05,dsg.s1=.1,mc.nsim=100,mc.nt=100,mc.pk=c(.7,.2,.05,.05),b0grid=c(.5,.9,20),b1grid=c(.7,1.7,20),b2grid=c(0,.8,20),uhat.bsave=rep(1,3),dir) {
# Description: performs forward simulation
# Input:
#   - dta.N: batch size. Posterior moments are obtained after every dta.N patients have been recruited into each treatment arm.
#   - th.u: prior for success probability theta ~ Beta(th.u,th.v)
#   - th.v: prior for success probability theta ~ Beta(th.u,th.v)
#   - th.au: hyperprior th.u ~ Gamma(th.au,th.bu) (mean is th.au/th.bu). Ignored if dsg.hyper and dsg.pool set to FALSE
#   - th.bu: hyperprior th.u ~ Gamma(th.au,th.bu) (mean is th.au/th.bu). Ignored if dsg.hyper and dsg.pool set to FALSE
#   - th.av: hyperprior th.v ~ Gamma(th.av,th.bv) (mean is th.av/th.bv). Ignored if dsg.hyper and dsg.pool set to FALSE
#   - th.bv: hyperprior th.v ~ Gamma(th.av,th.bv) (mean is th.av/th.bv). Ignored if dsg.hyper and dsg.pool set to FALSE
#   - dsg.uth: indicates what utility function to use. 0: use posterior mean in utility evaluation; 1: use true theta; 2: drop phaseIII term from utility
#   - dsg.hyper: 0: when simulating, treat hyper-param u,v as fixed; 1: simulate u,v from hyper-priors
#   - dsg.pool: 1: at the time of analysis pool via empirical Bayes; 0: don't pool (i.e. ignore hyper-prior)
#   - dsg.logsd: if set to TRUE the decision boundaries are linear in log(SD), if set to FALSE boundaries are linear in SD
#   - dsg.patlimit: maximum number of patients available at any time point. The default -1 means there's no limit.
#   - dsg-kini: initial number of treatments
#   - dsg-c0: coefficients for utility function
#   - dsg-c1: coefficients for utility function
#   - dsg-c2: coefficient for utility function. Pay-off for significant phase III result.
#   - dsg-p0: value of theta under H0 for phase III (e.g. success probability for standard of care)
#   - dsg-a: type I error for phase III study
#   - dsg-b: type II error for phase III study
#   - dsg-s0: min sd to terminate (log scale)
#   - dsg-s1: b1&b2 are the value of the boundaries when sd=dsg-s1
#   - mc.nsim: number of forward simulations
#   - mc.nt: time horizon in forward simulation expressed in number of batches.
#   - mc.pk: vector indicating probabilities of 0,1,2... new treatments arising at any point in time
#   - b0grid: grid on b0, specified as initial term, final term and number of points
#   - b1grid: grid on b0, specified as initial term, final term and number of points
#   - b2grid: grid on b0, specified as initial term, final term and number of points
#   - uhat.bsave: index for b0grid,b1grid,b2grid. For corresponding (b0,b1,b2) value extra information will be saved (trajectories, posterior moments etc.)
#   - dir: directory where output files will be written to
# Output: the following files will be written to the directory specified by 'dir' (defaults to working directory)
#   - blah blah blah


if (!is.logical(dsg.hyper)) stop('dsg.hyper must be TRUE or FALSE')
if (!is.logical(dsg.pool)) stop('dsg.pool must be TRUE or FALSE')
if (!is.logical(dsg.logsd)) stop('dsg.logsd must be TRUE or FALSE')
if (sum(mc.pk)!=1) stop('Elements of mc.pk must add up to 1')

dirini <- getwd()
if (!missing(dir)) { setwd(dir) }

write(paste("dta-N",dta.N,"\n"),"init.dsg")
write(paste("th-u",th.u),"init.dsg",append=TRUE)
write(paste("th-v",th.v),"init.dsg",append=TRUE)
write(paste("th-au",th.au),"init.dsg",append=TRUE)
write(paste("th-bu",th.bu),"init.dsg",append=TRUE)
write(paste("th-av",th.av),"init.dsg",append=TRUE)
write(paste("th-bv",th.bv,"\n"),"init.dsg",append=TRUE)
write(paste("dsg-uth",dsg.uth),"init.dsg",append=TRUE)
write(paste("dsg-hyper",as.numeric(dsg.hyper)),"init.dsg",append=TRUE)
write(paste("dsg-pool",as.numeric(dsg.pool)),"init.dsg",append=TRUE)
write(paste("dsg-logsd",as.numeric(dsg.logsd)),"init.dsg",append=TRUE)
write(paste("dsg-patlimit",as.numeric(dsg.patlimit),"\n"),"init.dsg",append=TRUE)
write(paste("dsg-kini",dsg.kini),"init.dsg",append=TRUE)
write(paste("dsg-c0",dsg.c0),"init.dsg",append=TRUE)
write(paste("dsg-c1",dsg.c1),"init.dsg",append=TRUE)
write(paste("dsg-c2",dsg.c2),"init.dsg",append=TRUE)
write(paste("dsg-p0",dsg.p0),"init.dsg",append=TRUE)
write(paste("dsg-a",dsg.a),"init.dsg",append=TRUE)
write(paste("dsg-b",dsg.b),"init.dsg",append=TRUE)
write(paste("dsg-s0",dsg.s0),"init.dsg",append=TRUE)
write(paste("dsg-s1",dsg.s1,"\n"),"init.dsg",append=TRUE)
write(paste("mc-nsim",mc.nsim),"init.dsg",append=TRUE)
write(paste("mc-nt",mc.nt),"init.dsg",append=TRUE)
cat("mc-Npk",length(mc.pk)-1,"\n",file="init.dsg",append=TRUE)
cat("mc-pk",mc.pk,"\n \n",file="init.dsg",append=TRUE)
cat("uhat-b0grid",b0grid,"\n",file="init.dsg",append=TRUE)
cat("uhat-b1grid",b1grid,"\n",file="init.dsg",append=TRUE)
cat("uhat-b2grid",b2grid,"\n",file="init.dsg",append=TRUE)
cat("uhat-bsave",uhat.bsave-1,"\n",file="init.dsg",append=TRUE)     #in C vectors start at position 0

.C("cforwsim",PACKAGE="seqdesphII")

setwd(dirini)

}
