read.forwsim <- function(dir) {
# Input:
#   - dir: directory containing output files from the forward simulation
# Output: a list containing the following elements
#   - b0: vector of values of b0
#   - b1: vector of values of b1
#   - b2: vector of values of b2
#   - bgrid: matrix containing all combinations of b0, b1 and b2 values
#   - u: mean utility (1 entry for each row in bgrid). Set to minimum value for those b0,b1,b2 not evaluated in the simul
#   - nreco: mean number of recommended treatments (1 entry for each row in bgrid)
#   - ntoreco: mean number of patients / number of recommended treatments (1 entry for each row in bgrid)
#   - pB: mean posterior probability that future trial gives significant result (1 entry for each row in bgrid)
#   - n3: mean sample size for future trial (1 entry for each row in bgrid)
#   - dt: mean number of patients tested with each treatment (1 entry for each row in bgrid)
#   - err: mean type I & II error rates (1 entry for each row in bgrid)
#   - mse: estimated mean square error (1 entry for each row in bgrid)
#   - traj: some sample trajectories
#   - term: terminal decisions and posterior moments for sample trajectories

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

# Read grid on (b0,b1,b2)
b0 <- scan('b0.sim',quiet=TRUE); b1 <- scan('b1.sim',quiet=TRUE); b2 <- scan('b2.sim',quiet=TRUE)
nb0 <- length(b0); nb1 <- length(b1); nb2 <- length(b2)
bgrid <- expand.grid(b0,b1,b2); names(bgrid) <- c('b0','b1','b2')
bgrid <- bgrid[order(bgrid[,1],bgrid[,2],bgrid[,3]),]

# Read expected utilities
u <- scan('U.sim',quiet=TRUE)
validb <- (bgrid$b1>bgrid$b0 & bgrid$b0>bgrid$b2)
u[!validb] <- min(u)

# Read other output (note: output files can have character values for not valid grid points, so read them as char)
# (note: ntoreco can have Inf values that under Windows can be saved in a format not readable by R. Fix: read as char, turn into numeric and NAs indicate where Inf values are)
nreco <- ntoreco <- pB <- n3 <- dt <- mse <- rep(0,nrow(bgrid))
err <- matrix(0,nrow=nrow(bgrid),ncol=2)

nreco[validb] <- as.numeric(scan('nreco.sim',quiet=TRUE,what='character')[validb])
ntoreco[validb] <- suppressWarnings(as.numeric(scan('ntoreco.sim',quiet=TRUE,what='character')[validb])); ntoreco[is.na(ntoreco)] <- Inf
pB[validb] <- as.numeric(scan('pB.sim',quiet=TRUE,what='character')[validb])
n3[validb] <- as.numeric(scan('n3.sim',quiet=TRUE,what='character')[validb])
dt[validb] <- as.numeric(scan('dt.sim',quiet=TRUE,what='character')[validb])  
err[validb,] <- as.numeric(cbind(scan('errI.sim',quiet=TRUE,what='character'),scan('errII.sim',quiet=TRUE,what='character'))[validb])
colnames(err) <- c('errI','errII'); err <- data.frame(err)
mse[validb] <- as.numeric(scan('mse.sim',quiet=TRUE,what='character')[validb])

#Files with additional results may be empty if index indicated some (b0,b1,b2) out of the grid. Use 'try' to avoid error message
traj <- try(read.table('traj.dsg'),silent=TRUE); if (mode(traj)!='character') names(traj) <- c('t','trt','m','logs','theta')
term <- try(read.table('term.sim'),silent=TRUE); if (mode(term)!='character') names(term) <- c('b0','b1','b2','sim','trt','m','logs','theta','term','t0','t1')

setwd(dirini)
return(list(b0=b0,b1=b1,b2=b2,bgrid=bgrid,u=u,nreco=nreco,ntoreco=ntoreco,pB=pB,n3=n3,dt=dt,err=err,mse=mse,traj=traj,term=term))

}
