iqp2sample3 <- function(x,y,times, error=1e-8, alpha=1, 
      w0=rep(1/length(x),length(x)), p0=rep(1/length(y), length(y)) ){
####
#### This function computes the empirical likelihood/NPMLE for 2 samples
#### under several equality constraints. 
#### x, y, --- the 2 sample observations, required inputs.
#### output (Weight) will be in the same order as input. 
#### If use optional input w0, p0, they need to be in the same order of x, y.
#### times, --- a set of give times used to get constraints. The 
#### constraints are: 1-F(times) = [1-G(times)]^alpha . 
#### If alpha is different from 1 then they are Lehmann alternative.
#### the optional input w0, p0 do not have to satisfy above constraints
#### but they need to be a valid probability mass function.
#### We assume   min(x) < times < max(x) and min(y) < times < max(y). 
#### Last changed 8/25/2000

###########  acquire original vales

  n <- length(x)            # smaple size
  m <- length(y) 

####    mat0a <- diag( 1/c(w0^2, p0^2) )     # original first entry
    mat0a <- diag( c(w0, p0) )               # the factorized version
    mat0b <- 1/c(w0, p0)                 # original second entry

    temp <- c(rep(1,n),rep(0,m))
    A1 <- cbind(temp, 1-temp)  
    A21 <- outer(x, times, FUN="<=")
    A22 <- outer(y, times, FUN="<=")
    A2 <- rbind(A21, -A22) 
    mat0c <- cbind(A1, A2)     # the constrain matrix, third entry

    di <- c(w0, p0) %*% A2 
    mat0d <- c(0,0,-di)     # original constrain values, forth entry

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

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

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

     while(diff > error) {               
####        mata <- diag(1/wp^2)      # the new first entry
        mata <- diag( wp )            # the factorized version
        matb <- 1/wp              # the new second entry
                                 # third entry do not need to be changed
        value <- solve.QP(mata, matb, mat0c, matd, meq=2+length(di), 
                         factorized=TRUE)
        wp <- wp + value$solution   # update the weight function
        diff <- sum(value$solution^2)    # new error 
     }

####     ratio <- exp(sum(log(n*wp)) )
  loglik <- sum(log(wp))
 
  list(WeightX=wp[1:n], WeightY=wp[(n+1):(n+m)], Value=loglik) 

}