iqpHybrid <- function(times1,d1,R1,times2,d2,R2, points, error=1e-9, eta=1){
#### This function computes the empirical likelihood/NPMLE for 2 samples
#### under several equality constraints. 
#### times12, d12, R12 --- the 2 sample observations, required inputs.
#### The times12 should be ordered, output (Weight) will be in the same order. 
#### d/R must be < 1, if the last value of (d/R) is 1, then delete it!!!
#### points, --- a set of given t's used to get constraints. The
#### constraints are: 1-F(t's) = [1-G(t's)]^eta .
#### 
#### We assume   min(times12) < points < max(times12).
#### Last changed 6/25/2001
    d12 <- c(d1, d2)
    R12 <- c(R1, R2)
    w012 <- d12/R12       # the Nelson-Aalen, unconstrained max location
    u012 <- log(1-w012)
    mat0a <- diag( (exp(-u012/2)-exp(u012/2))/sqrt(d12) )  
                             # first entry, the factorized version
    mat0b <- R12-d12 - (d12*exp(u012))/(1-exp(u012))  # second entry
    
    A21 <- outer(times1, points, FUN="<=")
    A22 <- outer(times2, points, FUN="<=")
    mat0c <- A2 <- rbind(A21, - eta*A22)  # the constrain matrix, third entry

    dd <- u012 %*% A2
    mat0d <- (-dd)     # original constrain values, forth entry
    k <- length(dd)

   value0 <- solve.QP(mat0a, mat0b, mat0c, mat0d, meq=k, factorized=TRUE) 
                                  # first calculation
    wp <- u012 + value0$solution     # the weight vector

########## now iteration

     matd <- rep(0, k)   # the new constrain values, forth entry
     diff <- 10          # original valus of accuracy

     while(diff > error) {               
        mata <- diag( (exp(-wp/2)-exp(wp/2))/sqrt(d12) )  
                           # first entry, the factorized version
        matb <- R12-d12 - (d12 * exp(wp))/(1-exp(wp))
                                          # the new second entry
                                 # third entry do not need to be changed
        value <- solve.QP(mata, matb, mat0c, matd, meq=k, factorized=TRUE)
        wp <- wp + value$solution   # update the weight function
        diff <- sum(value$solution^2)    # new error 
     }
##########
 loglik2 <- sum(d12*log(1-exp(wp))+(R12-d12)*wp)
 loglik1 <- sum(d12*log(w012)+(R12-d12)*u012)
   n <- length(times1)
   m <- length(times2)
  list(WeightX=wp[1:n], WeightY=wp[(n+1):(n+m)], ChiSq=2*(loglik1-loglik2))
}