####################
##   2019/01/30   ##
## Tatsuki Koyama ##
####################
##
## 3+3 design by simulation

three3 <- function(p, startingDose=1, oneThirdOk=FALSE){
        sq <- startingDose
        x.sq <- NULL
        n <- rep(0, length(p))
        x <- rep(0, length(p))
    stopNow <- FALSE
        nextDose <- startingDose
    while( !stopNow ){
        v <- rbinom(1,3, p[nextDose])
            n[ nextDose ] <- n[nextDose] + 3
            x[ nextDose ] <- x[nextDose] + v
                ptox <- x/n
        upDown <- -1 * (v>1) + 1 * (v==0) ## -1, 0, or 1
        nextDose <- nextDose + upDown
            sq <- c(sq, nextDose)
            x.sq <- c(x.sq, v)
        ## Stop if ...
        ## There is no dose level
        stopNow <- stopNow | (nextDose==0) | nextDose > length(p)
        if(nextDose>0){
        ## Already n=6 on the next dose
        stopNow <- stopNow | (n[nextDose] == 6)
        ## Already x=2 or 3 on the next dose
          tooMany <- 2 + (oneThirdOk)
        stopNow <- stopNow | (x[nextDose] >= tooMany)
        }
    }
        sq <- rev(rev(sq)[-1])
    ## getAnswer
    ToxP <- x/n
        eps <- ifelse(oneThirdOk, 0.01, 0)
        okDose <- ToxP < (1-0.001+eps)/3
        finalAnswer <- rev(which(okDose))[1]
            ## This will be NA if there is no answer.
        if(is.na(finalAnswer)) finalAnswer <- 0
        if(finalAnswer>0){
         ## if n<6 on the finalAnswer then additional n=3 on that dose
            nAtfinal <- n[finalAnswer]
            if(nAtfinal==3){ 
                    v <- rbinom(1,3, p[finalAnswer])
                    n[ finalAnswer ] <- n[ finalAnswer ] + 3
                    x[ finalAnswer] <- x[ finalAnswer ] + v
                    sq <- c(sq, finalAnswer)
                    x.sq <- c(x.sq, v)
                    }
    }
    nx.sq <- data.frame(DOSE=sq, TOX=x.sq)
    nx <- data.frame(N=n, TOX=x, PROB=round(x/n,2))
        finalAnswer <- ifelse(is.na(finalAnswer), 0, finalAnswer)
    list(nx.sq, nx, finalAnswer)
}

# three3(p=c(1/9,1/6,1/3,2/3), startingDose=1, oneThirdOk=!FALSE)

many3 <- function(B, p, startingDose, oneThirdOk){
    Answer <- N <- X<- numeric(B)
    for(i in 1:B){
        out <- three3(p, startingDose, oneThirdOk)
        Answer[i] <- out[[3]]
        N[i] <- sum(out[[2]][,1])
        X[i] <- sum(out[[2]][,2])
    }
    list(MTD=Answer,N=N,X=X)
    }

summary33 <- function(mny3){
    ## input is output of many 3
    list(
        MTD=table(mny3[[1]]), N=table(mny3[[2]]), X=table(mny3[[3]])
        )
}

# s1 <- many3(1000, p=c(0.01,0.05,0.1,1/3), 1, FALSE)
# s2 <- many3(1000, p=c(0.05,1/3,0.5,0.6), 1, FALSE)
# s3 <- many3(1000, p=c(1/3,1/3,1/3,1/3), 1, FALSE)
# summary33(s1)
# summary33(s2)
# summary33(s3)