tdms=function(X1,X2,m,epsilon,num.permutations=Inf){
## Input:
# X1 - network 1
# X2 - network 2
# m - minimum cluster size parameter
# epsilon - threshold parameter
# num.permutations - the number of random permutations
## Output:
# p - p value
# sN - observed test statistics
# mF1 - number of clusters in network 1
# number of clusters in network 2
s1=PLSnet(X1)
s2=PLSnet(X2)
n1=ncol(X1)
n2=ncol(X2)
X=cbind(X1,X2)
p=nrow(X)
n=n1+n2
F1=get.modules(abs(s1),m,epsilon)
F2=get.modules(abs(s2),m,epsilon)
G0=(F1>0)|(F2>0)
sNc=0
for (g in which(G0))
if ((F1[g]>0)&(F2[g]>0))
sNc=sNc+(F1[g]>0)*(F2[g]>0)*length(intersect(which(F1==F1[g]),which(F2==F2[g])))/length(union(which(F1==F1[g]),which(F2==F2[g])))
if (sum(G0)==0)
sN=0
else
sN=1-sNc/sum(G0)
R=choose(n,n1)
if (num.permutations<R){
permutation.list=sample(1:n,n1)
for (i in 2:num.permutations)
permutation.list=rbind(permutation.list,sample(1:n,n1))
}
else
permutation.list=comb(n,n1)
num.perm=nrow(permutation.list)
perm.sN=rep(0,num.perm)
i=1
perm.sG0=rep(0,num.perm)
mF1=max(F1)
mF2=max(F2)
perm.mF1=rep(0,num.perm)
perm.mF2=rep(0,num.perm)
perm.sF1=rep(0,num.perm)
perm.sF2=rep(0,num.perm)
while (i<=num.perm){
i1=as.vector(permutation.list[i,])
perm.X1=X[,i1]
perm.X2=X[,-i1]
perm.s1=PLSnet(perm.X1)
perm.s2=PLSnet(perm.X2)
perm.F1=get.modules(abs(perm.s1),m,epsilon)
perm.F2=get.modules(abs(perm.s2),m,epsilon)
perm.G0=(perm.F1>0)|(perm.F2>0)
sNc=0
for (g in which(perm.G0))
if ((perm.F1[g]>0)&(perm.F2[g]>0))
sNc=sNc+(perm.F1[g]>0)*(perm.F2[g]>0)*length(intersect(which(perm.F1==perm.F1[g]),which(perm.F2==perm.F2[g])))/length(union(which(perm.F1==perm.F1[g]),which(perm.F2==perm.F2[g])))
if (sum(perm.G0)==0)
perm.sN[i]=0
else
perm.sN[i]=1-sNc/sum(perm.G0)
i=i+1
}
p.value=mean(perm.sN>=sN)
list(p=p.value,sN=sN,mF1=mF1,mF2=mF2)
}