BLANKET=function(ranks1,ranks2,include.a=F){
 n=length(ranks1) 
 temp.matrix=outer(ranks1,ranks2,FUN="==") 
 num.list1=outer(ranks1,rep(1,n),FUN="*") 
 num.list2=outer(rep(1,n),ranks1,FUN="*") 
 temp1=apply(temp.matrix,1,cumsum) 
 a=apply(temp1,1,cumsum) 
 b=num.list1-a 
 c=num.list2-a 
 d=n-a-b-c
 exp.a=(num.list1/n)*num.list2
 exp.b=num.list1-exp.a
 exp.c=num.list2-exp.a
 exp.d=n-exp.a-exp.b-exp.c
 z=sqrt((a-exp.a)^2/(exp.a)+(b-exp.b)^2/(exp.b)+
      	(c-exp.c)^2/(exp.c)+(d-exp.d)^2/(exp.d))*sign(a-exp.a)
 pstat=1-pnorm(z)
 pfisher=1-phyper(a-1, num.list2, n-num.list2, num.list1)
 if (include.a==T){
      output=list(pstat,pfisher,a)
      names(output)=c("pstat","pfisher","a")
        }
 else {
      output=list(pstat,pfisher)
      names(output)=c("pstat","pfisher")
        }
 return(output)
 }
####### run an example
n=100
ranks1=1:n
ranks2=ranks1
lowerset=1:(n*.1)
ranks2[lowerset]=sample(ranks2[lowerset])
ranks2[-lowerset]=sample(ranks2[-lowerset])

date()
result=BLANKET(ranks1,ranks2)
date()
image(log10(result$pfisher))

### get empirical distribution of minima
num.perms=1000
num.common=c(400,500)
pmin.matrix=matrix(1,num.perms,length(num.common))
pmin20.matrix=matrix(1,num.perms,length(num.common))
pmin.a10.matrix=matrix(1,num.perms,length(num.common))

for (k in (1:length(num.common))){

pmin=rep(0,num.perms)
pmin20=rep(0,num.perms)
pmin.a10=rep(0,num.perms)

for (i in (1:num.perms)){
 ranks1=c(1:num.common[k])
 ranks2=sample(ranks1)
 result=BLANKET(ranks1,ranks2,include.a=T)
 pmin[i]=min(result$pfisher)
 pmin20[i]=min(result$pfisher[1:20,1:20])
 pmin.a10[i]=min(result$pfisher[result$a<=10])

 print(c(i,k))
 }

pmin.matrix[,k]=pmin
pmin20.matrix[,k]=pmin20
pmin.a10.matrix[,k]=pmin.a10

}

y1=apply(pmin.matrix,2,quantile,.05)
y2=apply(pmin20.matrix,2,quantile,.05)
y3=apply(pmin.a10.matrix,2,quantile,.05)

plot(num.common,y1,type="b",ylim=c(0,max(y1,y3,y2)))
lines(num.common,y2)
lines(num.common,y3,lty=3)


plot(num.common,apply(pmin20.matrix,2,quantile,.05))

plot(num.common,apply(pmin.a10.matrix,2,quantile,.05))

cbind(num.common,apply(pmin.matrix,2,quantile,.05),apply(pmin.a10.matrix,2,quantile,.05))




temp=read.table("C:\\Documents and Settings\\Administrator\\Desktop\\for_FWright_test_case.txt",sep="\t")
p1=temp[,2]
p2=temp[,4]
ranks1=rank(p1)
ranks2=rank(p2)
plot(ranks1,ranks2)

result=BLANKET(ranks1,ranks2)


###########

temp1=read.table("pmin.matrix")
temp2=read.table("pmin20.matrix")
temp3=read.table("pmin.a10.matrix")

temp1=temp1[,1:3]
temp2=temp2[,1:3]
temp3=temp3[,1:3]


threshold.all=c(apply(temp1,2,quantile, .05),apply(pmin.matrix,2,quantile,.05))
threshold.20=c(apply(temp2,2,quantile, .05),apply(pmin20.matrix,2,quantile,.05))
threshold.a10=c(apply(temp3,2,quantile, .05),apply(pmin.a10.matrix,2,quantile,.05))

plot(c(100,200,300,400,500),threshold.all,type="b",ylim=c(0,max(threshold.20)))
lines(c(100,200,300,400,500),threshold.20)
lines(c(100,200,300,400,500),threshold.a10,lty=2)

cbind(c(100,200,300,400,500),threshold.all,threshold.20,threshold.a10)



plot(c(100,200,300,400,500),log10(threshold.all),type="b")
lines(c(100,200,300,400,500),log10(threshold.20))
lines(c(100,200,300,400,500),log10(threshold.a10),lty=2)



****

#examples

BB=read.table("C:\\Users\\Kim\\Desktop\\BB_lists.txt",sep="\t")
BS=read.table("C:\\Users\\Kim\\Desktop\\BS_lists.txt",sep="\t")
Hammond=read.table("C:\\Users\\Kim\\Desktop\\Hammond_lists.txt",sep="\t")

par(mfrow=c(2,2))
BB.n=dim(BB)[1]
BB.ranks1=c(1:BB.n)
BB.ranks2=c(1:BB.n)[order(BB[,2])]
BB.ranks2=BB.ranks2[BB[,1]]
#cbind(BB.ranks1,BB.ranks2)[1:20,]
plot(BB.ranks1,BB.ranks2)

BS.n=dim(BS)[1]
BS.ranks1=c(1:BS.n)
BS.ranks2=c(1:BS.n)[order(BS[,2])]
BS.ranks2=BS.ranks2[BS[,1]]
plot(BS.ranks1,BS.ranks2)

Hammond.n=dim(Hammond)[1]
Hammond.ranks1=c(1:Hammond.n)
Hammond.ranks2=c(1:Hammond.n)[order(Hammond[,2])]
Hammond.ranks2=Hammond.ranks2[Hammond[,1]]
plot(Hammond.ranks1,Hammond.ranks2)

######

par(mfrow=c(2,2))
BB.result=BLANKET(BB.ranks1,BB.ranks2,include.a=T)
BB.p=c(min(BB.result$pfisher),
       min(BB.result$pfisher[BB.result$a<=10]),
       min(BB.result$pfisher[1:20,1:20]))
image(log10(BB.result$pfisher))

BS.result=BLANKET(BS.ranks1,BS.ranks2,include.a=T)
min(BS.result$pfisher[BS.result$a<=10])
BS.p=c(min(BS.result$pfisher),
       min(BS.result$pfisher[BS.result$a<=10]),
       min(BS.result$pfisher[1:20,1:20]))
image(log10(BS.result$pfisher))

Hammond.result=BLANKET(Hammond.ranks1,Hammond.ranks2,include.a=T)
min(Hammond.result$pfisher[Hammond.result$a<=10])
Hammond.p=c(min(Hammond.result$pfisher),
       min(Hammond.result$pfisher[Hammond.result$a<=10]),
       min(Hammond.result$pfisher[1:20,1:20]))
image(log10(Hammond.result$pfisher))


pvalue.matrix=rbind(BB.p,BS.p,Hammond.p)
dimnames(pvalue.matrix)[[2]]=c("min.all","min.a<=10","min.20X20")