### sum of the sum of log of rank statistics (SSLR)
SSLR=function(K) {
    SSLR1=0
    n=length(K)
    K=log(K)
    for (i in 1:n){
        SSLR1=SSLR1+sum(K[1:i])
    }
    SSLR1
}


### null distribution of SSLR
SSLRNull=function(n,N, no.sample=1000000){
    output=rep(0, no.sample)
    for (i in 1:no.sample) {
            output[i]= SSLR(sample(1:N, n))
    }
    output

}


#### permutation test of comparing multiple scoring functions, where Z is a matrix of the ranks by the scoring functions. 
compMdock=function(Z, no.perm=10000) {
    n=dim(Z)[[2]]
    output=matrix(NA, ncol=6, nrow=choose(n,2) )
    colnames(output)=c("function 1", "function 2","sslr 1", "sslr 2", "difference", "p value")
    count=1
    cn=colnames(Z)
    for (i in 1:(n-1)) {
        for (j in (i+1):n) {
            output[count, 1]=cn[i]; output[count, 2]=cn[j]
            temp=comp2dock(Z[,i], Z[,j], no.perm=no.perm)
            output[count, 3]=SSLR(Z[,i])
            output[count, 4]=SSLR(Z[,j])
            output[count, 5]=temp$dif
            output[count, 6]=temp$pvalue                         
            count=count+1
        }
    
    }
    output
}


### comparing two scoring functions, X and Y are the ranks of active compounds from two scoring functions. 
comp2dock=function(X,Y, no.perm=10000) {
    stat=SSLR(X)-SSLR(Y)
    stat.perm1=rep(0, no.perm)
    stat.perm2=rep(0, no.perm)
    KS=cbind(X, Y)
    for (i in 1:no.perm){
        B=rbinom(dim(KS)[[1]], 1,0.5)
        K1=B*KS[,1]+(1-B)*KS[,2]
        K2=B*KS[,2]+(1-B)*KS[,1]
        stat.perm1[i]=SSLR(K1)
        stat.perm2[i]=SSLR(K2)
    }
    p.value= mean(abs(stat.perm1-stat.perm2)>abs(stat))
    list(pvalue=p.value, dif=stat)
}


### rank sum test corresponds to AU-ROC
RS=function(K){
    RS1=0
    n=length(K)
    RS1=sum(K)
}


### null distribution of rank sum test
RSNull=function(n,N, no.sample=1000000){
    output=rep(0, no.sample)
    for (i in 1:no.sample) {
            output[i]= RS(sample(1:N, n))
    }
    output
}


####################
###Simulation Study   
####################

N=1000; n=10                ### N=total number of compounds;  n=number of active compounds
no.sample=1000000

### the null distribution of SSLR

SSLRNull1=SSLRNull(n,N,no.sample)
par(mar=c(5,5,5,5))
hist(SSLRNull1,nclass=200 , freq=F, xlab="SSLR Test", col="grey", main="", cex.lab=2, cex.axis=2)


### the following code compares SSLR with Rank Sum Test


### this function simulates the ranks of the active compounds
### "lambda" controls how early they are ranked and large "lambda" preferrs early recognition
ranksim=function(n,N,lambda){
    X=rep(0,n)
    while(length(unique(X))<n) {
        X=-log(1-sample(1:N, n)/N*(1-exp(-lambda)))/lambda*N
        X=round(X+0.5)
    }
    
    X
}


Truth=c(rep(1,n), rep(0, N-n))
no.iter=10000
lambda=c(0.0000000001, seq(0.5, 8, by=0.5))
rs=NULL                 ### ranks sum test 
sslr.u=NULL             ### sslr for unsorted ranks    
sslr.p=NULL             ### sslr for partially sorted ranks
sslr.i=NULL             ### sslr for completely sorted ranks

for (j in lambda){
    rs.test=rep(0, no.iter) 
    sslr.u.test=rep(0, no.iter)
    sslr.p.test=rep(0, no.iter)
    sslr.i.test=rep(0, no.iter)
    print(j)
    for (i in 1:no.iter) {
        X=sample(sort(ranksim(n,N,lambda=j)), prob=n:1) 
        Y=sample(setdiff(1:N, X))
        temp=rank(c(X,Y))[which(Truth==1)]
        rs.test[i]=RS(temp)
        sslr.u.test[i]=SSLR(sample(X))     ### ranks are unsorted
        sslr.p.test[i]=SSLR(X)             ### ranks are partially sorted 
        sslr.i.test[i]=SSLR(sort(X))       ### ranks are completely sorted
   
    }
    rs=cbind(rs, rs.test); sslr.u=cbind(sslr.u, sslr.u.test); sslr.p=cbind(sslr.p, sslr.p.test);sslr.i=cbind(sslr.i, sslr.i.test)
}


RSNull1=RSNull(n,N,no.sample)           ### null distribution of rank sum test
thresh.rs=quantile(RSNull1, 0.05)       ### threshold corresponds to alpha=0.05 type I error rate
thresh.sslr=quantile(SSLRNull1, 0.05)

result=cbind(colSums( rs<thresh.rs)/no.iter,colSums( sslr.u<thresh.sslr)/no.iter,colSums( sslr.p<thresh.sslr)/no.iter,colSums( sslr.i<thresh.sslr)/no.iter )
colnames(result)=c("RS", "SSLR.u", "SSLR.p", "SSLR.i")
rownames(result)=lambda

par(mar=c(5,5,2,2))
plot(lambda[1:17], result[1:17,1], col=1, pch = 20,type="o", xlab=expression(lambda), ylab="Power", cex.lab=1.2, cex.axis=1.2,cex=1.2)
points(lambda[1:17], result[1:17,2], col=2, pch = 20, type="o", lty=2,cex=1.2)
points(lambda[1:17], result[1:17,3], col=3,pch = 20,  type="o", lty=2,cex=1.2)
points(lambda[1:17], result[1:17,4], col=4, pch = 20, type="o", lty=2,cex=1.2)
leg.txt <- c("ROC", "SSLR.U", "SSLR.P", "SSLR.C")
legend("bottomright", leg.txt, col = c(1,2,3,4), pch = 20,bg = "white",lty=c(1,2,2,2,2), cex=1.2,inset = 0.01)