spower <- function(rcontrol, rinterv, rcens, nc, ni,
					test=logrank, nsim=500, alpha=.05, pr=TRUE) {

crit <- qchisq(1-alpha, 1)
group <- c(rep(1,nc), rep(2,ni))
nexceed <- 0

for(i in 1:nsim) {
  if(pr && i %% 10 == 0) cat(i,'')
  yc <- rcontrol(nc)
  yi <- rinterv(ni)
  cens <- rcens(nc+ni)
  y <- c(yc, yi)
  S <- cbind(pmin(y,cens), 1*(y <= cens))
  nexceed <- nexceed + (test(S, group) > crit)
}
nexceed/nsim
}
  

Quantile2 <- function(scontrol, hratio, 
					  dropin=function(times)0, 
					  dropout=function(times)0,
					  m=7500, tmax, qtmax=.001, mplot=200, pr=TRUE, ...) {

## Solve for tmax such that scontrol(t)=qtmax
dlist <- list(...)
k <- length(dlist) && !is.null(dlist)
f    <- if(k) function(x, scontrol, qt, ...) scontrol(x, ...) - qt
else function(x, scontrol, qt) scontrol(x) - qt

if(missing(tmax)) {
  if(k) tmax <- uniroot(f, c(0,1e9), scontrol=scontrol, qt=qtmax, ...)$root
	else tmax <- uniroot(f, c(0,1e9), scontrol=scontrol, qt=qtmax)$root
}
if(pr) cat('\nInterval of time for evaluating functions:[0,',
		   format(tmax),']\n\n')

## Generate sequence of times to use in all approximations and sequence
## to use for plot method

times <- seq(0, tmax, length=m)
tim   <- seq(0, tmax, length=mplot)
tinc  <- times[2]

## Approximate hazard function for control group
sc <- scontrol(times, ...)    
hc <- diff(-logb(sc))
hc <- c(hc, hc[m-1])/tinc  ## to make length=m

## hazard function for intervention group
hr <- rep(hratio(times), length=m)
hi <- hc*hr

## hazard for control group with dropin
di  <- rep(dropin(times),length=m)
hc2 <- (1-di)*hc + di*hi

## hazard for intervention group with dropout
do  <- rep(dropout(times),length=m)
hi2 <- (1-do)*hi + do*hc

## survival for intervention group
si  <- exp(-tinc*cumsum(hi))

## Compute contaminated survival function for control and intervention
sc2 <- if(any(di>0))exp(-tinc*cumsum(hc2)) else sc
si2 <- exp(-tinc*cumsum(hi2))


## Store all functions evaluated at shorter times vector (tim), for
## plotting
asing <- if(.R.)function(x)x else as.single
sc.p  <- asing(approx(times, sc,  xout=tim)$y)
hc.p  <- asing(approx(times, hc,  xout=tim)$y)
sc2.p <- asing(approx(times, sc2, xout=tim)$y)
hc2.p <- asing(approx(times, hc2, xout=tim)$y)

si.p  <- asing(approx(times, si,  xout=tim)$y)
hi.p  <- asing(approx(times, hi,  xout=tim)$y)
si2.p <- asing(approx(times, si2, xout=tim)$y)
hi2.p <- asing(approx(times, hi2, xout=tim)$y)

dropin.p  <- asing(approx(times, di, xout=tim)$y)
dropout.p <- asing(approx(times, do, xout=tim)$y)
hratio.p  <- asing(approx(times, hr, xout=tim)$y)
hratio2.p <- hi2.p/hc2.p

tim       <- asing(tim)

plot.info <- list(
  "C Survival"                   =list(Time=tim,Survival=sc.p),
  "I Survival"                   =list(Time=tim,Survival=si.p),
  "C Survival w/Dropin"          =list(Time=tim,Survival=sc2.p),
  "I Survival w/Dropout"         =list(Time=tim,Survival=si2.p),
  "C Hazard"                     =list(Time=tim,Hazard=hc.p),
  "I Hazard"                     =list(Time=tim,Hazard=hi.p),
  "C Hazard w/Dropin"            =list(Time=tim,Hazard=hc2.p),
  "I Hazard w/Dropout"           =list(Time=tim,Hazard=hi2.p),
  "Dropin"                       =list(Time=tim,Probability=dropin.p),
  "Dropout"                      =list(Time=tim,Probability=dropout.p),
  "Hazard Ratio"                 =list(Time=tim,Ratio=hratio.p),
  "Hazard Ratio w/Dropin+Dropout"=list(Time=tim,Ratio=hratio2.p))

## Create S-Plus functions for computing random failure times for
## control and intervention subject to dropin, dropout, and hratio

r <- function(n, what=c('control','intervention'), 
			  times, csurvival, isurvival) {
  what <- match.arg(what)
  approx(if(what=='control')csurvival else isurvival, 
		 times, xout=runif(n), rule=2)$y
}
asing <- if(.R.)function(x)x else as.single
formals(r) <- list(n=integer(0),
                   what=c('control','intervention'),
                   times=asing(times), csurvival=asing(sc2),
                   isurvival=asing(si2))

structure(r, plot.info=plot.info, 
		  dropin=any(di>0), dropout=any(do>0),
		  class='Quantile2')
}

print.Quantile2 <- function(x, ...) {
  attributes(x) <- NULL
  print(x)
  invisible()
}

plot.Quantile2 <- function(x, 
						   what=c('survival','hazard','both','drop','hratio',
							 'all'), dropsep=FALSE,
						   lty=1:4, col=1, xlim, ylim=NULL, 
						   label.curves=NULL, ...) {
  what <- match.arg(what)
  pi <- attr(x, 'plot.info')
  if(missing(xlim)) xlim <- c(0,max(pi[[1]][[1]]))
  dropin  <- attr(x, 'dropin')
  dropout <- attr(x, 'dropout')
  i <- c(1,2,if(dropin)3,if(dropout)4)
  if(what %in% c('survival','both','all')) {
	if(dropsep && (dropin|dropout)) {
	  labcurve(pi[1:2], pl=TRUE, lty=lty, col=col, xlim=xlim, ylim=ylim,
			   opts=label.curves)
	  labcurve(pi[i[-(1:2)]], pl=TRUE, lty=lty, col=col, xlim=xlim, ylim=ylim,
			   opts=label.curves)
	} else
	labcurve(pi[i], pl=TRUE, lty=lty, col=col, xlim=xlim, ylim=ylim,
			 opts=label.curves)
  }

  if(what %in% c('hazard','both','all')) {
	if(dropsep && (dropin|dropout)) {
	  labcurve(pi[5:6], pl=TRUE, lty=lty, col=col, xlim=xlim, ylim=ylim,
			   opts=label.curves)
	  labcurve(pi[4+i[-(1:2)]], pl=TRUE, lty=lty, col=col, xlim=xlim, ylim=ylim,
			   opts=label.curves)
	} else
	labcurve(pi[4+i], pl=TRUE, lty=lty, col=col, xlim=xlim, ylim=ylim,
			 opts=label.curves)
  }

  if(what=='drop' || (what=='all' && (dropin | dropout))) {
	i <- c(if(dropin)9, if(dropout)10)
	if(length(i)==0) i <- 10
    labcurve(pi[i], pl=TRUE, lty=lty, col=col, xlim=xlim, ylim=ylim,
			 opts=label.curves)
  }

  if(what %in% c('hratio','all')) {
    i <- c(11, if(dropin|dropout) 12)
    labcurve(pi[i], pl=TRUE, lty=lty, col=col, xlim=xlim, ylim=ylim,
			 opts=label.curves)
  }

  invisible()
}

logrank <- function(S, group) {
  y     <- S[,1]
  event <- S[,2]
  i     <- order(-y)
  y     <- y[i]
  event <- event[i]
  group <- group[i]
  x     <- cbind(group==1, group==2, (group==1)*event, (group==2)*event)
  s     <- rowsumFast(x, y, FALSE)
  nr1 <- cumsum(s[,1])
  nr2 <- cumsum(s[,2])
  d1  <- s[,3]
  d2  <- s[,4]
  rd  <- d1+d2
  rs  <- nr1+nr2-rd
  n   <- nr1+nr2
  oecum <- d1 - rd*nr1/n
  vcum  <- rd * rs * nr1 * nr2 / n / n / (n-1)
  sum(oecum)^2 / sum(vcum,na.rm=TRUE)
}

Weibull2 <- function(times, surv) {

z1 <- -logb(surv[1])
z2 <- -logb(surv[2])
t1 <- times[1]
t2 <- times[2]
gamma <- logb(z2/z1)/logb(t2/t1)
alpha <- z1/(t1^gamma)

g <- function(times, alpha, gamma) {exp(-alpha*(times^gamma))}
formals(g) <- list(times=NULL, alpha=alpha, gamma=gamma)
g
}

#Function to fit a Gompertz survival distribution to two points
#The function is S(t) = exp[-(1/b)exp(a+bt)]
#Returns a list with components a and b, and a function for
#generating S(t) for a vector of times

Gompertz2 <- function(times, surv) {

z1 <- logb(-logb(surv[1]))
z2 <- logb(-logb(surv[2]))
t1 <- times[1]
t2 <- times[2]
b  <- (z2-z1)/(t2-t1)
a  <- z1 + logb(b)-b*t1

g <- function(times, a, b) {exp(-exp(a+b*times)/b)}
formals(g) <- list(times=NULL, a=a, b=b)
g
}

Lognorm2 <- function(times, surv) {

z1 <- qnorm(1-surv[1])
z2 <- qnorm(1-surv[2])
sigma <- logb(times[2]/times[1])/(z2-z1)
mu    <- logb(times[1]) - sigma*z1

g <- function(times, mu, sigma) {1 - pnorm((logb(times)-mu)/sigma)}
formals(g) <- list(times=NULL, mu=mu, sigma=sigma)
g
}