## this test is a nonparametric test which extends friedman.test() to cases that allow replicates in blocks
## exact p value is calculated

TwoSampleBlockRank.test=function(y,groups, blocks, alternative="two.sided") {
   
    ## calculate statsitics
    
    stat.rank=0
    sizes=NULL
    ## remove blocks with missing cells
    index=NULL
    B=unique(blocks)
    G=unique(groups)
    for (i in B){
        if(length(unique(groups[blocks==i]))==1) index=c(index, i)
    } 
    
    
    if(!is.null(index)) {
        for (i in index) {blocks=ifelse(blocks==i, NA, blocks)}
        temp=!is.na(blocks)
        blocks=blocks[temp]
        groups=groups[temp]
        y=y[temp]
    }
    
    
    B=unique(blocks)
    G=unique(groups)
    
    lb=length(B)
    for (i in 1:lb){
       stat.rank= sum(rank(y[blocks==B[i]]) [groups[blocks==B[i]]==G[1]])+stat.rank
       sizes=c(sizes, choose(sum(blocks==B[i]), sum(groups[blocks==B[i]]==G[1])))
    }
    ## NULL Distribution
    ranksum.matrix=matrix(NA, ncol=max(sizes), nrow=lb)
    
    for (i in 1:length(B)) {
        temp=rowSums(choose.enumerate(rank(y[blocks==B[i]]), n=sum(groups[blocks==B[i]]==G[1])))
        ranksum.matrix[i,1:length(temp) ]= temp
    }
    
    b.length=NULL
    for(i in 1:dim(ranksum.matrix)[[1]]) {
       b.length=c(b.length, sum(!is.na(ranksum.matrix[i, ])))    
    }
    temp=sort(b.length, index.return=T)
    min.length=min(b.length)
    b.length=b.length[temp$ix]
    ranksum.matrix=matrix(ranksum.matrix[temp$ix,], nrow=lb)
  
    ranksum=NULL
    continue1=T
    index=rep(1, length(B))
    while(continue1){
        ranksum=c(ranksum, colSums(matrix(ranksum.matrix[, 1:min.length], ncol=min.length,nrow=lb)))
        if(lb==1) continue1=F
        else if(index[lb]==b.length[lb]) {
             continue2=T
             temp=lb-1
             while(continue2) {
                 if(temp==1) {continue2=F; continue1=F}
                 else if(index[temp]==b.length[temp]) {continue2=T; temp=temp-1}
                 else{
                        continue2=F; index[temp]=index[temp]+1; a=ranksum.matrix[temp,1];
                        ranksum.matrix[temp, 1:(b.length[temp]-1)]=ranksum.matrix[temp, 2:(b.length[temp])]
                        ranksum.matrix[temp, b.length[temp]]=a;
                        index[(temp+1):length(index)]=1
                 }
                
            }
        }
        else{
            index[lb]=index[lb]+1
            a=ranksum.matrix[lb,1];
            ranksum.matrix[lb, 1:(b.length[lb]-1)]=ranksum.matrix[lb, 2:(b.length[lb])]
            ranksum.matrix[lb, b.length[lb]]=a        
        }
    }
    
    if(alternative=="two.sided") {
        p.value=2*mean(ranksum>=stat.rank) 
        if(p.value>1)    p.value=2*mean(ranksum<=stat.rank)    
        if(p.value>1)    p.value=1                             
    }
    if(alternative=="greater") p.value=mean(ranksum>=stat.rank)
    if(alternative=="less") p.value=mean(stat.rank>=ranksum)
    output=list(p.value=p.value, stat=stat.rank)

}


choose.enumerate=function(x, n=NA){
    
    if(is.na(n)) stop ("The number to choose is not specified!")    
    lx=length(x)
    nrow=choose(lx, n)
    output=matrix(0, ncol=n, nrow=nrow)
    continue1= T
    count1=1
    index=1:n
    while(continue1) {
        output[count1, ]=x[index]
        count1=count1+1; 
        if(count1>nrow) continue1=F
        else {
            if(index[n]==length(x)) {
                continue2=T
                temp=n-1
                max.limit=lx-1
                while(continue2){
                    if(index[temp]==max.limit) {continue2=T; temp=temp-1; max.limit=max.limit-1}
                    else{continue2=F;index[temp]=index[temp]+1; index[temp:n]=seq(index[temp],index[temp]+(n-temp) )  }
                }
            }
            else {index[n]=index[n]+1}
         }
    }
    output
}