```{r,echo=FALSE,message=FALSE,results=FALSE}
opts_chunk$set(echo=FALSE)
options(warn=-1)
```

```{r}
source("functionsRB.R")

library(isotone)
probFun <- function(dose) pnorm((dose-0.729)/0.084)
#	function(x) 1/(1+exp(-(-15.36 + 21.02*x)))
#probFun <- function(x) 1/(1+exp(-1*x))
#pts <- seq(from=0,to=1,by=.1)
#plot(pts,probFun(pts))

estED <- function(x, qtile,maxTurns,probFun) {
	nobs <- nrow(x)
	res <- rep(NA,maxTurns)
  nturns <- x[nobs,2]
	kp <- x[,2] <= maxTurns
#	if (sum(kp) > maxN) { fail <- TRUE; break }
 	isofit <- isoreg(x[kp,-2])
  ordx <-  isofit$x
  fits <- isofit$yf	
# 	print(summary(cbind(ordx,fits)))
    if (length(unique(fits)) == 1) last(x[,1]) else approx(isofit$yf[kp], ordx[kp], xout=qtile,rule=2)$y
 }

# BCD (ED = .90)

target <- .5
#target <- .9
targetDose <- qnorm(target)*.084 + .729 
if (round(10*target)==5) {
startDoses <- c(.4,.5,.6,.7)
maxN <- 30
} else {
startDoses <- c(.2,.5,.8,1)
maxN <- 60
}

nTurns <- c(6,10)
#target <- .90
downProb <- 1/target - 1
chgVec <- c(1,0,-1)
mcN <- 1000
stepSizes <- c(.025,.05,.075,.1,.125,.15,.25)
```


```{r}
ssize <- integer(mcN) 
#sampleSizes <- array(dim=c(mcN,length(startDoses),length(stepSizes)),
#	   dimnames= list(NULL, as.character(startDoses),as.character(stepSizes)))
samples <- list()
doseRanges <- list()
NAs <- rep(NA,maxN)
#starts <- array(dim=c(length(startDoses),length(stepSizes)),
	    #dimnames=list(start=as.character(startDoses),step=as.character(stepSizes)))
set.seed(1010101)
for (startDose in startDoses) {
  for (stepSize in stepSizes) {
#startDose <- startDoses[2]
#stepSize <- stepSizes[3]
# startDose <- startDoses[1] ; stepSize <- stepSizes[1] 
doses <- startDose + stepSize * seq(from=-10,to=100,by=1)
doses <- doses[doses > 0]
startLevel  <- which(doses == startDose)
doseRanges[[paste("Start:",startDose," Step:", stepSize)]] <- doses
nDoses <- length(doses)
probResp <- probFun(doses)
unifArr <- array(runif(maxN*mcN),dim=c(maxN,2,mcN))
unifArr[,2,] <-  unifArr[,2,] <  downProb
result <- array(dim=c(maxN,3,mcN))
for (run in 1:mcN) {
  unifVec <- unifArr[,1,run]
  downVec <- -1 * unifArr[,2,run]
  prevDose <- curDose <- startLevel
  reversals <- 0
  yVec <-  doseVec <- revVec <- NAs
  lastMove <- 0
  lastY <- NA
  for (i in 1:maxN) {
      doseVec[i] <- curDose
      yVec[i] <- yi <- ( unifVec[i] < probResp[curDose] )
#      if (reversals == stopReverse) break
      nextStep <- switch(as.character(yi),"FALSE"=1,"TRUE"= downVec[i])
      nextDose <- min(max(1,curDose + nextStep),nDoses)
      revVec[i] <- reversals <- reversals +  (!is.na(lastY) & (lastY != yi))
      curDose <- nextDose
      lastY <- yi
  }
  result[,,run] <-  cbind(doses[doseVec],revVec,yVec)
}
samples[[paste("Start:",startDose," Step:", stepSize)]] <-  result
}
}
```


```{r}

estimates <- list()
for (nTurni in nTurns) {
estimates[[as.character(nTurni)]] <- lapply(samples, function(x) { apply(x,3,estED,target, maxTurns = nTurni) } )
}

dimn <-  c(list(what=names(unisum(rnorm(100)))), 
	 lapply(list(step=stepSizes, start=startDoses, nTurn = nTurns), as.character))
sumryEst <- array(unlist(lapply(estimates, lapply,unisum),recursive=TRUE),
	  dim=sapply(dimn,length),dimnames=dimn)
emses	 <- apply(sumryEst,2:4,function(x) sqrt((x[1]-targetDose)^2 + x[2]^2))  
	    
sumryProb <-   array(unlist(lapply(lapply(estimates,lapply,probFun),lapply,unisum),
	  recursive=TRUE), dim=sapply(dimn,length),dimnames=dimn)
pmses	 <- apply(sumryProb,2:4,function(x) sqrt((x[1]-target)^2 + x[2]^2))  
#rscpmses <- pmses*(sqrt(2*pi)*.084)/exp(-(qnorm(target)^2)/2)
    
ssNeeded <- function(samples,nTurns,maxN) {
sszs <- list()
for (resi in seq(along=samples)) {
ss <- apply( samples[[resi]], 3, function(x) { min(which(x[,2] >= nTurns))  } ) 
ss[is.infinite(ss)] <- maxN
sszs[[resi]]  <- ss
}
sszs
}

sszs <- list()
for (nTurni in nTurns) sszs[[as.character(nTurni)]] <- ssNeeded(samples,nTurns=nTurni,
	   maxN=maxN)

dimn <-  lapply(list(rep=1:mcN,step=stepSizes, start=startDoses, nTurns = nTurns), as.character)

ssAry <- array(unlist(sszs,recursive=TRUE),
	  dim=sapply(dimn,length),dimnames=dimn)
ssDf <- array.to.data.frame(ssAry,data.name="ssize")
save(estimates,samples,emses,pmses,ssDf, file=paste("upDownSim.rda",sep=""))

```

```{r}

for (starti in dimn[[3]]) {
	   dfi <- ssDf[ssDf$start==starti,]
	   bxp <- with(dfi,{ par(oma=c(0,1,3,0)) ;
	   	twoway.boxplot(ssize,nTurns,step, oma=c(0,0,9,0)) } )
	   mtext(line=.5,at=bxp[[1]],round(t(emses[,starti,]),3),cex=.8)
	   mtext(line=.5,at=-.3,"RMSE(Est):",cex=.9,adj=1)
	   mtext(line=1.5,at=bxp[[1]],round(t(pmses[,starti,]),3),cex=.8)
	   mtext(line=1.5,at=-.3,"RMSE(Prob):",cex=.9,adj=1)
	   mtext(line=3,paste("Starting Dose =",starti),cex=1.2)
	   mtext(line=1,side=1,at=-.3,"Turns",cex=.9,adj=1)
	   mtext(line=2,side=1,at=-.3,"Step Size",cex=.9,adj=1)
}

```