discemlik <- function(x, d, K, fun, tola=.Machine$double.eps^.25, theta)
{
#  The function discemlik() calculates the jumps of H(t) and the log empirical
#  likelihood ratio for right censored data under the constraint  
#  sum_i[g(ti, theta)log(1- dH(ti))] = K  where H(t) is the cumulative hazard. 
#  The likelihood function is the discrete version, see Cox and Oakes book.
#  When K is a const, this is an implicit function for theta. 
#  But when theta held fix and K change this becomes an explicit function 
#  for K.  Ref: Zhou (1999) U of Kentucky Tech report 36?. 
#  Or see Fang, Hui(2000) P.26.
#  Discrete version, for implicit/explicit functions. One sample.
#  
#  works with Splus version 3.4 and R version 1.0 and later + survival5. 
#  It uses the function survfit() so make sure you have them.
#  Written by Mai Zhou (mai@ms.uky.edu)  Feb. 1999 version 1.0
#
#  INPUTS:
#     x        Observed (possibly censored) times:  x=min(t,C).
#     d        Indicator of censoring: d=1 not censor, d=0 censor.
#     K        the constant that the integrated hazard is = to under H0. 
#              For mediam, fun g(t, theta)=I[t<theta] and K=log2. 
#              For explicit functions, theta considered fixed and K changes. 
#     fun      Fuction of 2 variables, g(t, theta), to compute 
#              the constraint. It has to be able to take a vector input 
#              in the first variable and a scalar (parameter) in the second.
#     tola     Optional tolerance value in solving the lambda = 0 equation.
#     theta    Value of the implicit parameter constraint, passed to fun g.
#
#  OUTPUTS:
#     logemlik   -2Log empirical likelihood ratio (ALR or Poisson version)
#     lambda     Lagrange multiplier
#     times      the location of the jumps, could missing last one
#     jumps      the size of jumps, the last jump=1 may be missing
#     nits       Number of iterations used
#
#   Examples of the input function fun is:  (indicator function) 
#                 fun2 <- function(x, theta) { as.numeric(x <= theta) }
# 
#   bug:  theta value has to be in certain range for this to work properly.
#         Assumes fun() is non-negative.
#   example:  x <- c(1, 2, 3, 4, 5, 6, 5, 4, 3, 4, 1, 2.4, 4.5) 
#             d <- c(1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1)       then 
# > discemlik(x,d,K=-0.7,fun2,theta=4)    # (fun2 as above) should return 
#                                           a list including (among others)
#      $"discrete.-2logemlikRatio"   
#       [1] 0.1446350
#########################################################################
n <- length(x) 
if(n <= 2) stop("Need more observations")
if(length(d) != n ) stop("length of x and d must agree") 
if(any((d!=0)&(d!=1))) stop("d must be 0/1's for censor/not-censor")
if(!is.numeric(x)) stop("x must be numeric values --- observed times")
#############################################################
#x <- pmin(xx,dd)                #THis is for simulations. Delete the 2 lines
#d <- as.numeric(xx <= dd)       #and change xx to x , dd to d  

temp <- summary(survfit(Surv(x,d),se.fit=F,type="fleming",conf.type="none"))
otime <- temp$time         # only uncensored time?  Yes. 
orisk <- temp$n.risk 
odti <- temp$n.event

###if the last jump is of size 1, we need to drop last jump from computation
last <- length(orisk)    
if (orisk[last] == odti[last]) {
                     otime <- otime[-last] 
                     orisk <- orisk[-last]
                     odti  <- odti[-last]
                     }
######## compute the function g(ti, theta) 
gti <- fun(otime,theta) 

### the constrain function. To be solved in equation later.

constr <- function(x, Konst, gti, rti, dti, n) { 
                  rtiLgti <- rti + x*n*gti
                  sum(gti*log((rtiLgti - dti)/rtiLgti) ) -  Konst } 

##############################################################

differ <- constr(0, Konst=K, gti=gti, rti=orisk, dti=odti, n=n)

if( abs(differ) < tola ) { lam <- 0 } else {
    step <- 0.5/sqrt(n) 
    if(abs(differ) > 200*log(n)*step )   #Why 60 ? 
    stop("given theta value is too far away from theta0")

    mini<-0
    maxi<-0            
######### assume the constrain function is increasing in lam (=x) 
    if(differ > 0) { 
    mini <- -step 
    while(constr(mini, Konst=K, gti=gti, rti=orisk, dti=odti, n=n) > 0 &&
                  mini > -200*log(n)*step ) 
    mini <- mini - step 
    } 
    else { 
    maxi <- step 
    while(constr(maxi, Konst=K, gti=gti, rti=orisk, dti=odti, n=n) < 0 &&
                  maxi < 200*log(n)*step ) 
    maxi <- maxi+step 
    } 

    if(constr(mini, Konst=K, gti=gti, rti=orisk, dti=odti, n=n)*constr(maxi, 
                Konst=K, gti=gti, rti=orisk, dti=odti, n=n) > 0 )
    stop("given theta/K is/are too far away from theta0/K0")

# Now we solve the equation, to get lambda, to satisfy the constraint of Ho

    temp2 <- uniroot(constr,c(mini,maxi), tol = tola, 
                  Konst=K, gti=gti, rti=orisk, dti=odti, n=n)  
    lam <- temp2$root 
} 
##########################################################################
rPlgti <- orisk + n*lam*gti

loglik <- 2*sum(odti*log(rPlgti/orisk) + 
             (orisk-odti)*log(((orisk-odti)*rPlgti)/(orisk*(rPlgti-odti)) ) ) 

#?is that right? YES the -2log lik ratio. 
# Notice the output time and jumps has less the last point.
list("discrete.-2logemlikRatio"=loglik, lambda=lam, times=otime, 
                jumps=odti/rPlgti)
}