power=function(N,a,f,pi,theta=NA,t0,type="CSH",alpha=0.05,
KMev0=NA,KMcr0=NA,KMev1=NA,KMcr1=NA,
CIFev0=NA,CIFcr0=NA,CIFev1=NA,CIFcr1=NA)
{
##################################################################
#### CALCULATES POWER
## for difference between two groups
## need to import the function eqsolve
## n=total number of patients which will be accrued
## theta=hazard ratio:
##   if type=CSH then is the ratio of the cause specific hazrard
##   if type=HCIF then is the ratio of teh hazards of CIF
## a=accrual time, 
## f=follow-up time, 
## t0=time at which KM and CIF are given
## a,f and t0 are all in the same units
## CIFev0=CIF for event in standard group at t0
## CIFcr0=CIF for competing risks in standard group at t0
## CIFev1=CIF for event in experimental group at t0
## CIFcr1=CIF for competing risks in experimental group at t0
## KMev0=KM estimate for free of event in standard group at t0
## KMcr0=KM estimate for free of event in standard group at t0
## KMev1=KM estimate for free of event in experimental group at t0
## KMcr1=KM estimate for free of event in experimental group at t0
## pi=proportion in standard group, default=0.5
## alpha= alpha level, default=0.05
##### ASSUMPTIONS #####
## exponential distribution for both groups and both events
## independence between the types of events
## the competing risks the same in the standard and competing risks
####################################################################
if (!is.na(KMev0) & !is.na(KMcr0))
  {
lev0=-log(KMev0)/t0
lcr0=-log(KMcr0)/t0
  }
if (is.na(KMev0) & is.na(KMcr0) & !is.na(CIFev0) & !is.na(CIFcr0))
  {
CIF0=CIFev0+CIFcr0
part=-log(1-CIF0)/(t0*CIF0)
lev0=CIFev0*part
lcr0=CIFcr0*part
  }

if (is.na(lev0) | is.na(lcr0)) stop("Not enough information given")
lev1=NA
lcr1=lcr0

if (!is.na(KMev1) & !is.na(KMcr1))
  {
lev1=-log(KMev1)/t0
lcr1=-log(KMcr1)/t0
  }
if (is.na(KMev1) & is.na(KMcr1) & !is.na(CIFev1) & !is.na(CIFcr1))
  {
CIF1=CIFev1+CIFcr1
part=-log(1-CIF1)/(t0*CIF1)
lev1=CIFev1*part
lcr1=CIFcr1*part
  }

if (is.na(lev1) & is.na(theta)) stop("theta is missing")
if (type=="CSH" & !is.na(theta)) lev1=lev0/theta
if (type=="HCIF" & !is.na(theta)) CIFev1=1-(1-CIFev0)**(1/theta)
if (is.na(lev1)) lev1=eqsolve(CIFev1,lcr1,t0)
l0=lev0+lcr0
l1=lev1+lcr1

if (is.na(theta) & type=="CSH") theta=lev0/lev1
if (is.na(CIFev0)) CIFev0=lev0/l0*(1-exp(-l0*t0))
if (is.na(CIFev1)) CIFev1=lev1/l1*(1-exp(-l1*t0))
if (is.na(theta) & type=="HCIF") theta=log(1-CIFev0)/log(1-CIFev1)
pev0=lev0/l0*(1-(exp(-l0*f)-exp(-l0*(a+f)))/(l0*a))
pev1=lev1/l1*(1-(exp(-l1*f)-exp(-l1*(a+f)))/(l1*a))
pev=pev0*pi+pev1*(1-pi)
nev=N*pev

zalpha=qnorm(1-alpha/2)
zbeta=sqrt(nev*pi*(1-pi))*abs(log(theta))-zalpha
power=pnorm(zbeta)
return(power)
}

eqsolve=function(fev,lcr,t0)
{
lev1=0.0001
lev2=2
k=0
f1=fev-lev1/(lcr+lev1)*(1-exp(-(lcr+lev1)*t0))
f2=fev-lev2/(lcr+lev2)*(1-exp(-(lcr+lev2)*t0))
if (f1*f2>0) stop("Either lambda is much too high or CIF is much too small")
fm=10
while (abs(fm)>0.000001)
  {
if (k>10000)
    {
stop("No convergence")
    }
levm=lev1+(lev2-lev1)/2
fm=fev-levm/(lcr+levm)*(1-exp(-(lcr+levm)*t0))
if (f1*fm<0) 
     {
lev2=levm
f2=fm
k=k+1
     }
if (f2*fm<0) 
     {
lev1=levm
f1=fm
k=k+1
     }
  }

r=round(levm,4)
return(r)
}