optbound <- function(forwsim,errImax,errIImax,util='u',max=TRUE,avedraws=TRUE,smooth=FALSE,span=.75) {
# Description: finds optimal boundaries
# Input:
#   - forwsim: list containing forward simulation results as returned by read.forwsim
#   - errImax: upper bound of type I error probability. Only boundaries satisfying it will be included in the optimization
#   - errIImax: upper bound of type II error probability.
#   - util: name of the element in forwsim containing the expected utility that is to be optimized
#   - max: if TRUE the goal is to maximize the utility, otherwise it is to minimize it
#   - avedraws: when set to TRUE, multiple optimal solutions will be averaged to give a single solution.
#   - smooth: should the utility surface be loess smoothed previous to optimization?
#   - span: the parameter alpha which controls the degree of smoothing as defined in 'loess'
# Ouput: list containing the following elements
#   - b: values of b0,b1,b2 maximizing the expected utility
#   - u: estimated expected utility at the optimum
#   - err: matrix with estimated type I & II error probabilities at the optimum in 1st and 2nd columns, respectively
#   - dt: mean number of patients tested with each treatment at the optimum
#   - mse: estimated mean square error in parameter estimation at the optimum
#   - span: if smooth set to TRUE, value of the smoothing parameter 

u <- forwsim[[util]]
validb <- (forwsim$bgrid$b1>forwsim$bgrid$b0 & forwsim$bgrid$b0>forwsim$bgrid$b2)

#Smoothing
if (smooth) {
  if (length(span)>1) {
    r2 <- rep(NA,length(span)); maxr2 <- 0
    test <- rbinom(length(u),1,p=1/3)
    test[forwsim$bgrid$b0==min(forwsim$b0) | forwsim$bgrid$b0==max(forwsim$b0) | forwsim$bgrid$b1==min(forwsim$b1) | forwsim$bgrid$b1==max(forwsim$b1) | forwsim$bgrid$b2==min(forwsim$b2) | forwsim$bgrid$b2==max(forwsim$bgrid$b2)] <- 0
    u.test <- data.frame(forwsim$bgrid[test==1,],u=u[test==1]); u.learn <- data.frame(forwsim$bgrid[test==0,],u=u[test==0])
    for (i in 1:length(span)) {
       u.smooth <- loess(u~b0+b1+b2,data=u.learn,control=loess.control(surface='interpolate',statistics='approximate',trace.hat='approximate'),span=span[i],degree=1)
       u.smooth <- predict(u.smooth,u.test[,1:3],se=FALSE)
       r2[i] <- 1 - sum((u.test$u - u.smooth)^2)/sum((u.test$u-mean(u.test$u))^2)
     }
    span <- span[r2==max(r2)]
  }
  u[validb] <- loess(u[validb]~forwsim$bgrid[validb,1]+forwsim$bgrid[validb,2]+forwsim$bgrid[validb,3],control=loess.control(surface='interpolate',statistics='approximate',trace.hat='approximate'),span=span,degree=1)$fitted
  u[!validb] <- min(u,na.rm=TRUE)
} else { span <- NA }

#Find optimal b value
if (min(forwsim$err[validb,1])<=errImax) {                                   #if some rules satisfy typeI error constraint
  s <- validb & (forwsim$err[,1]<=errImax)
  b <- forwsim$bgrid[s,]; u <- u[s]; err <- forwsim$err[s,]; dt <- forwsim$dt[s]; mse <- forwsim$mse[s]
  if (min(err[,2])<=errIImax) {                                              #if some rules satisfy typeII error constraint
    s <- (err[,2]<=errIImax)
    b <- b[s,]; u <- u[s]; err <- err[s,]; dt <- dt[s]; mse <- mse[s]
    bopt <- if (max) { (u==max(u)) } else { (u==min(u)) }
  } else { bopt <- (err[,2]==min(err[,2])) }                                 #if no rules satisfy typeII error constraint
} else {  bopt <- validb & (forwsim$err[,1]==min(forwsim$err[,1])) }         #if no rules satisfy typeI error constraint
 
b <- b[bopt,]; u <- u[bopt]; err <- err[bopt,]; dt <- dt[bopt]; mse <- sqrt(mse[bopt])
if ((nrow(b)>1) & avedraws) { b <- colMeans(b); u <- mean(u); err <- colMeans(err); dt <- mean(dt); mse <- mean(mse) }

return(list(b=as.numeric(b),u=u,err=err,dt=dt,mse=mse,span=span))

}
