iqp2sample <- function(x,y,times, error=1e-8,
      w0=rep(1/length(x),length(x)), p0=rep(1/length(y), length(y)) ){

#### x, y, times  ------  are required input.
#### 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.
#### constraints are: F(times)==G(times). 
#### 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). 
#### Also, there must be some x, and y in between each times. 

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

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

    mat0a <- diag( 1/c(w0^2, p0^2) )     # original first entry
    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)) 
                                  # first calculation
    wp <- c(w0, p0) + value0$solution     # the weight vector

########## 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
        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))
        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) 

}