##### 	Program below generates a single trial that uses design from 
	Ivanova " Dose-finding for continuous and ordinal outcomes 
	with a monotone objective function: A unified approach". 
	Scenarios are from Cheung and Chapell (Biometrics, 2000).
	The program is written for binary outcomes but can be easily modified to use other distribution. 


nonpart=function(xy)
{
		# input data 
model=xy[1]	# dose-response scenario 
ga=xy[2]	# target quantile 
n=xy[3]		# total sample size
eps=xy[4]	# desgin parameter delta
grsize=xy[7]	# cohort size
startup=xy[8]	# if start-up is desired, startup=1; otherwise no start-up

startupsize=round(log(.5)/log(1-ga)+.1,0)


if (model==1)	{truetox=c(.05,.10,.20,.30,.50,.70)		#pr(tox)
		}
if (model==2)	{truetox=c(.30,.4,.52,.61,.76,.87)		#pr(tox)
		}
if (model==3)	{truetox=c(.05,.06,.08,.11,.19,.34)		#pr(tox)
		}
if (model==4)	{truetox=c(.06,.08,.12,.18,.40,.71)		#pr(tox)
		}
if (model==5)	{truetox=c(.00,.0,.03,.05,.11,.22)		#pr(tox)
		}

le=length(truetox)
trials=rep(0,le)
y=rep(0,le)   		# toxicity outcome, y=1 if toxicity, y=0 if not 
dose0=1			#starting dose 
dosenum=dose0
count=0

if (startup==1) {while (sum(y)==0 && sum(trials)<startupsize*(le-dose0+1) ) {	trials[dosenum]=trials[dosenum]+startupsize
										resp=rbinom(1,startupsize,truetox[dosenum])
										y[dosenum]=y[dosenum]+resp
										dosenum=ifelse(resp>0,max(dosenum-1,1),min(le,dosenum+1))
				  		     				}
	        }

count=sum(trials)
while (count<n)	{groupsize1=min(n-count,grsize)
		 trials[dosenum]=trials[dosenum]+groupsize
		 resp=rbinom(1,groupsize,truetox[dosenum])
		 y[dosenum]=y[dosenum]+resp
		 qhat=y[dosenum]/trials[dosenum]
		 qhat1=max(qhat,0.001)
		 tstat=(qhat-ga)/sqrt(qhat1*(1-qhat1))*sqrt(trials[dosenum])
		 if (trials[dosenum]==1) dosenum1=dosenum
		 if (trials[dosenum] >1) {	dosenum1=dosenum
						if (tstat< (-eps)) dosenum1=dosenum+1
						if (tstat> (eps)) dosenum1=dosenum-1
		 			 }
		 dosenum=dosenum1
		 dosenum=max(dosenum,1)
		 dosenum=min(dosenum,le)
		 count=count+groupsize
		}
## estimating the maximum tolerated dose
trind=ifelse(trials>0,1,0)*(1:le)
if (trials[1]==0) {trind=trind[-1]; dob=1} else dob=0
qstar=isot(c(y[trind],trials[trind]))
####  modify qstar so that the target dose is estimated according for the algorithm from Ivanova 
for (j in 1:length(qstar)) if (qstar[j]<ga) qstar[j]=qstar[j]+1/1000/3^(7-j)
ebestd=order(round(abs(qstar-ga),5))[1] + dob
ebestd.v=rep(0,le)
ebestd.v[ebestd]=1
c(ebestd.v,trials,y,model,ga,n,eps)
}

nonpart(c(1,.25,24,.1,.1,1,1,1,3))