Split <- function(data,by){
    tmp <- as.list(seq_along(data))
    names(tmp) <- names(data)
    by <- eval(substitute(by),tmp,parent.frame())
    if(inherits(by,"formula")){
        if(length(by)>2){
            warning("Extranous lhs in 'by' argument")
            by <- by[-2]
        }
        by <- all.vars(by)
    }
    by <- data[by]
    j <- multique(by)
    split.data.frame(data,j)
}

complete_csets <- function(formula,
                           data,
                           weights,
                           idname="id",
                           choicename="choice"){
    altv <- deparse(formula[[2]])
    has_weights <- FALSE
    if(!missing(weights)){
        weights <- deparse(substitute(weights))
        has_weights <- TRUE
        }
    by <- formula[-2]
    dby <- data[all.vars(by)]
    j <- multique(dby)
    data[["__group__"]] <- j
    splt <- Split(data, by = by)
    if (has_weights) 
        splt <- lapply(splt, complete_csets1w, altv, weights)
    else
        splt <- lapply(splt, complete_csets1, altv, choicename)
    r <- lapply(splt,nrow)
    res <- do.call(rbind,splt)
    i <- res[["__set__"]]
    j <- res[["__group__"]]
    ij <- multique(i,j)
    res[["__set__"]] <- NULL
    res[["__group__"]] <- NULL
    res[[idname]] <- ij
    res
}


complete_csets1 <- function(data, altv, choicename){
    nms <- names(data)
    a <- data[[altv]]
    data[[altv]] <- NULL
    i <- muq(a)
    ua <- attr(i,"unique")
    n <- length(ua)
    j <- multique(data)
    r <- length(attr(j,"unique"))
    mat <- matrix(0,nrow=n,ncol=r)
    mat[cbind(i,j)] <- 1
    d <- data[!duplicated(j),,drop=FALSE]
    i <- as.vector(row(mat))
    j <- as.vector(col(mat))
    res <- d[j,,drop=FALSE]
    res[[choicename]] <- as.vector(mat)
    res[[altv]] <- ua[i]
    nms <- c(nms,choicename)
    res <- res[nms]
    res[["__set__"]] <- j
    res
}

complete_csets1w <- function(data, altv, weights) {
    nms <- names(data)
    a <- data[[altv]]
    w <- data[[weights]]
    data[[altv]] <- NULL
    data[[weights]] <- NULL
    i <- muq(a)
    ua <- attr(i, "unique")
    n <- length(ua)
    j <- multique(data)
    r <- length(attr(j, "unique"))
    mat <- matrix(0, nrow = n, ncol = r)
    mat[cbind(i, j)] <- w
    d <- data[!duplicated(j), , drop = FALSE]
    i <- as.vector(row(mat))
    j <- as.vector(col(mat))
    res <- d[j, , drop = FALSE]
    res[[weights]] <- as.vector(mat)
    res[[altv]] <- ua[i]
    res <- res[nms]
    res[["__set__"]] <- j
    res
}


muq <- function(x,sort=FALSE,drop.na=TRUE){
    u <- unique(x)
    if(drop.na)
        u <- na.omit(u)
    if(sort)
        u <- sort(u)
    structure(match(x,u),
              unique=u)
}


multique <- function(x,...){
    if(is.list(x))
        data <- x
    else
        data <- list(x,...)
    m <- length(data)
    j <- muq(data[[1]],drop.na=FALSE)
    if(m > 1){
        r <- length(attr(j,"unique"))
        for(k in 2:m){
            j_k <- muq(data[[k]],drop.na=FALSE)
            j <- j + r*(j_k-1)
            j <- muq(j)
            r <- length(attr(j,"unique"))
        }
    }
    j
}
