######################################
# Function to simulate survival data #
# n = sample size		     #
# lambda = hazard rate 		     #
# atime = accrual time (months)      #
# nsim = number of data sets	     #
######################################


simarray <- function(n, lambda, atime, nsim)  {

	simdat <- array(NA, dim=c(n, 5, nsim)) 

		for (i in 1:nsim) {
		s.time <- round(rexp(n,lambda))
		accrual.time <- round(runif(n,0,atime))
		fu.time <- 60-accrual.time
		event <- ifelse(fu.time < s.time,0,1)
		s.time <- ifelse(fu.time < s.time, fu.time, s.time)
		id <- as.integer(rank(accrual.time, ties.method = "first"))
		mat <- cbind(id, s.time, accrual.time, fu.time, event)
		simdat[,,i] <- mat[order(accrual.time),]
		}

	assign("simdat", simdat, env=.GlobalEnv)

}

###########################################
# Function to calculate the MLE of lambda #
# lmba: sequence of values for lambda	  #
###########################################

getlik <- function(lmba, simdat) {

	likvals <- matrix(NA, ncol=dim(simdat)[3], nrow=length(lmba))
 
	likmat <- matrix(NA, ncol=dim(simdat)[3], nrow = length(lmba))
	mle <- likatmle <- rep(NA, dim(simdat)[3])

	for (i in 1:dim(simdat)[3]) {
		
		keep <- ifelse(is.na(simdat[,1,i]),0,1)
		data <- simdat[keep==1,,i]
		mle[i] <- 1/(sum(data[,5]*(data[,2]) + (1 - data[,5])*data[,4]))*sum(data[,5])
		likatmle[i] <- round(sum(-1*(data[,5]*(log(1/mle[i]) + (data[,2])*mle[i]) + 
				(1 - data[,5])*mle[i]*data[,4])), 2)

		for (j in 1:length(lmba)) {
			likmat[j,i] <- round(sum(-1*(data[,5]*(log(1/lmba[j]) + (data[,2])*lmba[j]) + 
			(1 - data[,5])*lmba[j]*data[,4])), 2)
					}

		likvals[,i] <- likmat[,i]
		likmat[,i] <- likmat[,i] - likatmle[i]
		}

	#assign("likmat", likmat, env=.GlobalEnv)
	#assign("mle", mle, env=.GlobalEnv)

	return(list("likmat"=likmat,"mle"=mle))

}

################################################################
# Get 1/8 support intervals - point where lik function = 0.125 #
################################################################

getints <- function(likmat, mle, lmba, cutoff=8,time=NULL) {

low <- upp <- rep(NA, length(mle))

for (i in 1:length(mle)) {

low[i] <- ifelse( is.na(likmat[1,i]), 0, min(lmba[exp(likmat[,i]) >= 1/cutoff]))
upp[i] <- ifelse( is.na(likmat[1,i]), 1, max(lmba[exp(likmat[,i]) >= 1/cutoff]))


}

likpars <- cbind(mle, low, upp)

if(is.null(time)==F) likpars <- exp(-likpars*time)

assign("likpars", likpars, env=.GlobalEnv)

}


###################################################
# Format data for the number of people to look at #
# Uses data through end of follow-up		  #
# num: Number of patients to evaluate		  #
###################################################

getlookdat <- function(simdat, num, cuttime) {

lookdat <- array(NA, dim=c(num, dim(simdat)[2], dim(simdat)[3]))

for (i in 1:dim(simdat)[3]) {
timdiff <- simdat[,4,i][simdat[,1,i]<=num] - min(simdat[,4,i][simdat[,1,i]<=num])
fu.timen <- timdiff + cuttime
eventn <- ifelse(fu.timen < simdat[,2,i][simdat[,1,i]<=num],0,1)
s.timen <- ifelse(fu.timen < simdat[,2,i][simdat[,1,i]<=num], fu.timen, simdat[,2,i][simdat[,1,i]<=num])

lookdat[,,i] <- cbind(simdat[,1,i][simdat[,1,i]<=num], s.timen, simdat[,3,i][simdat[,1,i]<=num], fu.timen, eventn)
}

return(lookdat)

}

##############################################################
# Plot mle and support intervals for each simulated data set #
# The first function (plotmle) plots data at the end of trial#
# The 2nd function (plotbounds) plots data at each look	     #
# all.looks: array of the data for every 5 pts. - see below  #
##############################################################

plotmle <- function(likpars, color, lsim) {
plot(range(likpars), range(0:dim(likpars)[1]), type="n", xlab="Survival Probability at 6 Months", ylab="Simulated Data Set")
points(likpars[,1], 1:dim(likpars)[1], pch = 19)
segments(likpars[,2], 1:dim(likpars)[1], likpars[,3], 1:dim(likpars)[1], col=color)
abline(v = lsim, lty=2, col="grey")
}

plotbounds <- function(all.looks, color, lsim, xlab) {
par(bg="black", fg="white", col.axis = "white", col.lab="white")
plot(range(all.looks[,2:3,]), range(0:dim(all.looks)[1]), type="n", xlab=xlab, ylab="Simulated Data Set")
for (i in 1:dim(all.looks)[3]) {
segments(all.looks[,2,i], 1:dim(all.looks)[1], all.looks[,3,i], 1:dim(all.looks)[1], col=color[i])
}
abline(v = lsim, lty=2, col="grey")
}



#####################################################################
# Function to simulate survival data from a mixture of two exponentials #
# n = sample size		     #
# lambda = hazard rate 		     #
# atime = accrual time (months)      #
# nsim = number of data sets	     #
######################################


mixture.exp <- function(n1, lambda1, n2, lambda2, atime, nsim, study.time)  {

      n <- n1+n2
	simdat <- array(NA, dim=c(n, 5, nsim)) 
		for (i in 1:nsim) {
		s.time1 <- round(rexp(n1,lambda1))
		s.time2 <- round(rexp(n2,lambda2))
		s.time <- c(s.time1, s.time2)
		accrual.time <- round(runif(n,0,atime))
		fu.time <- study.time-accrual.time
		event <- ifelse(fu.time < s.time,0,1)
		s.time <- ifelse(fu.time < s.time, fu.time, s.time)
		id <- as.integer(rank(accrual.time, ties.method = "first"))
		mat <- cbind(id, s.time, accrual.time, fu.time, event)
		simdat[,,i] <- mat[order(accrual.time),]
		}
	return(simdat)
}

###########################################
# cut dataset after every kth person
# reaches 'cuttime' months
###########################################

interimlooks<-function(simdat, k, cuttime) {

lookdat <- array(NA, dim=c(dim(simdat)[1], dim(simdat)[2], dim(simdat)[3]))

for (i in 1:dim(simdat)[3]) {
	id <- simdat[,1,i]
	atime <- simdat[,3,i]
	ftime <- simdat[,4,i]
	etime <- simdat[,2,i]
	event <- simdat[,5,i]

	time <- atime[k]+cuttime
	keep <- ifelse(atime<=time,1,0)
	etime.new <- apply(cbind(etime, time - atime),1,min) 
	event.new <- ifelse(atime + etime < time, event, 0)
	ftime <- time - atime	

	n <- sum(keep)

	lookdat[1:n,,i] <- cbind(id, etime.new, atime, ftime, event.new)[keep==1,]
	}

return(lookdat)

}

###############
# make general function for getting mle and interval
#############
lik.result<-function(lmba, data, cutoff=8, time=NULL) {

liks<-getlik(lmba,data)
likpars <- getints(liks$likmat, liks$mle, lmba, time=time)
#likpars[,2][likpars[,1]==0] <- 0 # Replaces interval to cover whole line in case of no events #
#likpars[,3][likpars[,1]==0] <- 1
return(likpars)
}

#############
# function for getting x month survival time and 95% CI
#############

freq.int<-function(data,time=time) {
	st<-Surv(data[,1],data[,2]) 
	fit<-survfit(st)
	est<-min((summary(fit)$surv)[summary(fit)$time<=time])
	lower<-min((summary(fit)$lower)[summary(fit)$time<=time])
	upper<-min((summary(fit)$upper)[summary(fit)$time<=time])
	return(c(est,lower,upper))
}

###################
# plot lik and freq results after all pts.
###################

plot.both <- function(lik,freq,true) {

	xl<-range(lik,freq)
	yl<-c(-5,nrow(lik))
	plot(xl,yl,type="n",xlab="% surviving",ylab="dataset index")
	for(i in 1:nrow(lik)) {
		lines(lik[i,c(2,3)],rep(i-0.1,2),col=2)
		points(lik[i,1],i-0.1,col=2,pch=16)
		lines(freq[i,c(2,3)],rep(i+0.1,2),col=4)
		points(freq[i,1],i+0.1,col=4,pch=16)
	}
	abline(v=true)
	abline(h=0.5,lty=3)
	lines(range( lik[,1]),rep(-4, 2),col=2,lwd=2)
	lines(range(freq[,1]),rep(-2, 2),col=4,lwd=2)
title("Red = Likelihood; Blue = K-M Estimate")

return(NULL)
}

######################
# simulate "true" 1 year survival prob rate from mixture of exponentials
#####################

simtrue<-function(N,h1,h2,time, study.time) {

	est<-0
	data<-mixture.exp(1000,h1,1000,h2,1, N, study.time)
	for(i in 1:N) {
		sd<-Surv(data[,2,i], data[,5,i])
		fit<-survfit(sd)
		est<-est+min((summary(fit)$surv)[summary(fit)$time<=time])/N
	}
	
	return(est)
}

##############################
# generate data with assumed % cured
##############################

mixture.cured <- function(pct, lambda, n, atime, nsim, study.time)  {

	n1<-round(n*(1-pct))
	n2<-n-n1
	simdat <- array(NA, dim=c(n, 5, nsim)) 
		for (i in 1:nsim) {
		accrual.time <- round(runif(n,0,atime))
		fu.time <- study.time-accrual.time
		s.time1 <- round(rexp(n1,lambda))
		s.time <- c(s.time1, fu.time[(n1+1):n]+1)
		event <- ifelse(fu.time < s.time,0,1)
		s.time <- ifelse(fu.time < s.time, fu.time, s.time)
		id <- as.integer(rank(accrual.time, ties.method = "first"))
		mat <- cbind(id, s.time, accrual.time, fu.time, event)
		simdat[,,i] <- mat[order(accrual.time),]
		}
	return(simdat)
}

######################
# simulate "true" 1 year survival prob rate from mixture of exponential and cured
########################

simtrue.cured<-function(pct,N,h1, study.time) {

	est<-0
	data<-mixture.cured(pct, h1, 1000, 1, N, study.time)
	for(i in 1:N) {
		sd<-Surv(data[,2,i], data[,5,i])
		fit<-survfit(sd)
		est<-est+min((summary(fit)$surv)[summary(fit)$time<=time])/N
	}
	
	return(est)
}