.packageName <- "RankAggreg"
`BruteAggreg` <-
function(x, k, weighted=FALSE, index.weights=NULL, distance=c("Spearman", "Kendall")){
    distance <- match.arg(distance, c("Spearman", "Kendall"))
    x <- x[,1:k]
    distinct <- apply(x, 1, function(y) ifelse(length(unique(y)) < k, 1, 0))
    if (sum(distinct) >= 1)
        stop("Elements of Each Row Must Be Unique")
    if (nrow(x)<2)
        stop("X must have more than 1 row") 
    
    if(weighted){
        index.weights <- index.weights[,1:k]
        #standardize weights:
        index.weights <- t(apply(index.weights,1,function(z) (z-min(z))/(max(z)-min(z))))
        if(index.weights[1,k]!=0)
            index.weights <- 1-index.weights
        if(dim(x)[1] != dim(index.weights)[1] || dim(x)[2] != dim(index.weights)[2])
        stop("Dimensions of x and weight matrices have to be the same")
    }

       
    comp.list <- unique(sort(as.vector(x)))
    n <- length(comp.list)

    if (k > n)
        stop("k must be smaller or equal to n") 

    library(gtools)
    perms <- permutations(n,k,comp.list)
    if(distance=="Spearman")
        f.y <- spearman(x, perms, index.weights, weighted)
    else
        f.y <- kendall(x, perms, index.weights, weighted)
    list(top.list=perms[which.min(f.y),],fn.score=min(f.y), distance=distance)
}

`RankAggreg` <-
function(x, k, index.weights=NULL, use.weights=FALSE, method=c("CrossEntropy", "GeneticAlgorithm"), distance=c("Spearman", "Kendall"), 
            rho=.01, weight=.5, N=5*k*length(unique(sort(as.vector(x)))), error = .001, maxIter = 100,
            popSize=100, CP=1, MP=.001, informative=FALSE, v1=NULL, verbose=TRUE)
{
   
    method <- match.arg(method, c("CrossEntropy", "GeneticAlgorithm"))
    distance <- match.arg(distance, c("Spearman", "Kendall"))
    x <- x[,1:k]
    distinct <- apply(x, 1, function(y) ifelse(length(unique(y)) < k, 1, 0))
    if (sum(distinct) >= 1)
        stop("Elements of Each Row Must Be Unique")
    if (nrow(x)<2)
        stop("X must have more than 1 row") 

    compr.list <- unique(sort(as.vector(x)))
    n <- length(compr.list)
    
    if(method=="CrossEntropy"){
        comp.list <- 1:n
        x <- t(apply(x,1, function(xx) match(xx,compr.list)))}

    if(use.weights){
        index.weights <- index.weights[,1:k]
        #standardize weights:
        index.weights <- t(apply(index.weights,1,function(z) (z-min(z))/(max(z)-min(z))))
        if(index.weights[1,k]!=0)
            index.weights <- 1-index.weights
        if(dim(x)[1] != dim(index.weights)[1] || dim(x)[2] != dim(index.weights)[2])
        stop("Dimensions of x and weight matrices have to be the same")
    }

    if (k > n)
        stop("k must be smaller or equal to n") 

   
    if(method=="CrossEntropy"){
        v <- array(dim=c(n,k,maxIter))
        if(informative)
            v[,,1] <- as.matrix(v1)
        else
            v[,,1] <- 1/(n)
        y <- vector("numeric")
        diff <- vector("numeric")

        Nhat <- round(rho*N)
        if (Nhat < 5)
            stop("rho is too small")

        t <- 1
        repeat
        {
            cands <- mcmcProc(v[,,t], N, comp.list=comp.list)
            
            if(verbose)
                if(distance=="Spearman")
                    cat(" Calculating Spearman scores... \n")
                else
                    cat(" Calculating Kendall scores... \n")
                    
            if(use.weights)
                if(distance=="Spearman")
                    f.y <- spearman(x, cands, index.weights, weighted=TRUE)
                else
                    f.y <- kendall(x, cands, index.weights, weighted=TRUE)
            else
                if(distance=="Spearman")
                    f.y <- spearman(x, cands, weighted=FALSE)
                else
                    f.y <- kendall(x, cands, weighted=FALSE)
                    
            if(verbose)
                cat("Iteration", t, "Done! ")

            fy <- sort(f.y, ind=TRUE)
            y[t] <- fy$x[Nhat]
            good.cand <- cands[f.y <= y[t],]
            
            
            v[,,t+1] <- upd.prob(good.cand, v[,,t], weight, comp.list)
            
            best.cand <- compr.list[cands[fy$ix[1],]]
            rm(cands) # clean up
            
            y.l <- paste(best.cand, sep="", collapse=",")

            diff[t] <- sum(abs(v[,,t+1]-v[,,t]))/(n*k)
            
            if(verbose)     
                cat("\n",    c("Optimal value: ", y[t],
                        "\n Optimal List:  ",y.l, 
                       " \n Prob. Diff:    ",  diff[t]), 
                      ifelse(diff[t] < error,"<",">"), error, "\n\n")   

            if(diff[t] < error)          
                break 

            t <- t + 1
            if (t > 200){
            cat("Did not converge after 200 iterations. Please increase sample size N")
            break}
        }
    } else{
        #generate initial population randomly
        cands <- matrix(0, popSize, k)
        for(i in 1:popSize)
            cands[i,] <- sample(compr.list, k)
            
        #calculate obj. fn
        if(use.weights)
            if(distance=="Spearman")
                f.y <- spearman(x, cands, index.weights, weighted=TRUE)
            else
                f.y <- kendall(x, cands, index.weights, weighted=TRUE)
        else
            if(distance=="Spearman")
                f.y <- spearman(x, cands, weighted=FALSE)
            else
                f.y <- kendall(x, cands, weighted=FALSE)
        
        best.cand <- cands[which.min(f.y),]
        bestevery <- min(f.y)    
        
        conv=FALSE
        t <- 1
        iter <- 0
        while(!conv){
            #selection probability
            minf <- min(f.y)
            p.y <- (max(f.y)+1-f.y)/sum((max(f.y)+1-f.y))
            cpy <- cumsum(p.y)
            
            #select cands for the next generation
            ind <- runif(popSize)
            ind2 <- rep(0, popSize)
            for(i in 1:popSize)
                ind2[i] <- sum(ind[i] > cpy)+1
            cands <- cands[ind2,]

            # cross-over
            pairstocross <- floor(popSize*CP/2)
            samp <- sample(1:popSize, pairstocross*2)
            pointsofcross <- sample(2:k, pairstocross, replace=TRUE)
            for(i in 1:pairstocross){ 
                for(j in pointsofcross[i]:k){
                # this loop performs partially matched crossover (PMX) described in Section 10.5 of 
                # Data Mining: Concepts, Models, Methods, and Algorithms by Mehmed Kantardzic (2003)

                    t1 <- cands[samp[i],j]
                    t2 <- cands[samp[i+pairstocross],j]
                    
                    if(!is.na(t3 <- match(t2, cands[samp[i],])))
                        cands[samp[i], t3] <- t1
                    if(!is.na(t3 <- match(t1, cands[samp[i+pairstocross],])))
                        cands[samp[i+pairstocross], t3] <- t2
                                                
                    cands[samp[i], j] <- t2
                    cands[samp[i+pairstocross], j] <- t1
                }
            }              
                          
            # random mutations with probability MP
            mutations <- round(popSize*k*MP)
            rows <- sample(1:popSize, mutations, replace=TRUE)
            cols <- sample(1:k, mutations, replace=TRUE)
            for(i in 1:mutations)
                if(length(temp <- setdiff(compr.list, cands[rows[i],]))!=0)
                    cands[rows[i], cols[i]] <- sample(temp,1)
            
            
            #calculate obj. fn
            if(use.weights)
                if(distance=="Spearman")
                    f.y <- spearman(x, cands, index.weights, weighted=TRUE)
                else
                    f.y <- kendall(x, cands, index.weights, weighted=TRUE)
            else
                if(distance=="Spearman")
                    f.y <- spearman(x, cands, weighted=FALSE)
                else
                    f.y <- kendall(x, cands, weighted=FALSE)
  
            y.l <- paste(cands[which.min(f.y),], sep="", collapse=",")
            if(verbose)     
                cat("\n", "Iteration", t, ": ",  c("Optimal value: ", min(f.y),
                        "\n Optimal List:  ", y.l, "\n"))
            
            if(minf == min(f.y))
                iter <- iter+1
            else
                iter <- 1
            
            if(iter == 5 || t >= maxIter)
                conv=TRUE
            
            if(min(f.y) < bestevery){
                best.cand <- cands[which.min(f.y),]
                bestevery <- min(f.y)
            }
            t <- t+1
        }
    }
    list(top.list=best.cand, optimal.value=ifelse(method=="CrossEntropy", fy$x[1], bestevery), 
    sample.size = ifelse(method=="CrossEntropy", N, popSize), num.iter=t, method=method, distance=distance)
}

`kendall` <-
function(x, y, weights=NULL, weighted=TRUE)
{
    # Inputs:   x - lists to be combined
    #           y - candidate lists
    #           weights - weight matrix if weights to be used
    #           weighted - boolean if weights to be used
    # Outputs:  f.y - distances between each of candidate list and 
    #               x measured by Spearman formula
    
    k <- ncol(y)
    N <- nrow(x)

    # Kendall distances: regular and weighted
    kendall.dist <- function(x, y)
    {
        K=0
        n <- length(x)
        for (i in 1:(n-1))
            for (j in i:n)
                if((x[i] > x[j] & y[i] < y[j]) | (x[i] < x[j] & y[i] > y[j]))
                    K=K+1
        K
    }
    
    kendall.dist.w <- function(x, y, weight, k)
    {   
        K=0
        n <- length(x)
        for (i in 1:(n-1))
            for (j in (i+1):n)
                if((x[i] > x[j] & y[i] < y[j]) | (x[i] < x[j] & y[i] > y[j]))
                    K=K+abs(weight[ifelse(x[i]>k,k,x[i])]-weight[ifelse(y[j]>k,k,y[j])])
        K
    }


    if(weighted){   
        f.y <- apply(y, 1, function(z){
        sum <- 0
        for(i in 1:N){
            ul <- unique(c(z,x[i,]))
            rank.z <- match(ul,z,nomatch=k+1)
            rank.x <- match(ul,x[i,],nomatch=k+1)
            sum <- sum + kendall.dist.w(rank.x, rank.z, weights[i,], k)
        }
        sum
        })}
    else{
        f.y <- apply(y, 1, function(z){
            sum(apply(x,1,function(q){
                ul <- unique(c(z,q))
                rank.z <- match(ul,z,nomatch=k+1)
                rank.q <- match(ul,q,nomatch=k+1)
                kendall.dist(rank.z,rank.q)
            }))
        })
    }       
        
}

`mcmcProc` <-
function(v, N, thin=20, burn.in=round(N*thin/100), comp.list, verbose=TRUE)
{
    n <- nrow(v)
    k <- ncol(v)
    M <- (N+burn.in)*thin

    samples <- matrix(0,N+burn.in,k)
    ind <- 1
    last <- comp.list[sample(1:n,n)]    
    
    vs <- sample(1:k, M, replace=TRUE)
    hs <- sample(1:n, M, replace=TRUE)
    unifs <- runif(M)   
    
    if(verbose)
        cat("MCMC: ")   
    for (i in 1:M)
    {      
        x <- last
        vert <- vs[i]
        horz <- hs[i]

        temp <- last[vert]      
        x[vert] <- last[horz]
        x[horz] <- temp
    
        if(horz > k)        
            A <- v[x[vert],vert]/v[last[vert],vert]
        else
            A <- v[x[vert],vert]*v[x[horz],horz]/
                (v[last[vert],vert]*v[last[horz],horz])
    
        if (A > unifs[i]){
            if(i%%thin==0)      
                samples[i/thin,] <- x[1:k]
            last <- x
        }
        else
            if(i%%thin==0)
                samples[i/thin,] <- last[1:k]

        if(i%%(M/10)==0 && verbose)
            cat(paste(i*100/M, "% ",sep=""))
    }
    samples[-(1:burn.in),]
}

`spearman` <-
function(x, y, weights, weighted=TRUE)
{   
    # Inputs:   x - lists to be combined
    #           y - candidate lists
    #           weights - weight matrix if weights to be used
    #           weighted - boolean if weights to be used
    # Outputs:  f.y - distances between each of candidate list and 
    #               x measured by Spearman formula
    
    k <- ncol(y)
    N <- nrow(x)

    ### Spearman footrule functions (regular and weighted)
    subtr <- function(x,y)
        sum(abs(x-y))

    subtrw <- function(x,y,weight,k){
        # puts a weight of 0 for all x > k (weight[k]=0 due to normalizaton)
        sum(abs(weight[ifelse(x>k,k,x)]-weight[ifelse(y>k,k,y)])*abs(x-y))
        }
    ####
    
    if(weighted){   
        f.y <- apply(y, 1, function(z){
        sum <- 0
        for(i in 1:N){
            ul <- unique(c(z,x[i,]))
            rank.z <- match(ul,z,nomatch=k+1)
            rank.x <- match(ul,x[i,],nomatch=k+1)
            sum <- sum + subtrw(rank.x, rank.z, weights[i,], k)
        }
        sum
        })}
    else{
        f.y <- apply(y, 1, function(z){
            sum(apply(x,1,function(q){
                ul <- unique(c(z,q))
                rank.z <- match(ul,z,nomatch=k+1)
                rank.q <- match(ul,q,nomatch=k+1)
                subtr(rank.z,rank.q)
            }))
        })
    }       
        
}

`upd.prob` <-
function(samples, v, weight, comp.list)
{
    p <- matrix(0, nrow=nrow(v), ncol=ncol(v))
    s <- nrow(samples)
    p <- apply(samples,2,function(x) 
        table(x)[match(as.character(comp.list), dimnames(table(x))[[1]])]/s)
    p[is.na(p)] <- 0
    (1-weight)*v + weight*p
}

