tdcind=function(X1,X2,num.permutations=Inf){
## Input:
# X1 - network 1
# X2 - network 2
# num.permutations - the number of random permutations
## Output:
# p - p value
# d - observed test statistics
s1=plsnet(X1)
s2=plsnet(X2)
n1=ncol(X1)
n2=ncol(X2)
X=cbind(X1,X2)
n=n1+n2
nG=nrow(X)
di=rep(0,nG)
for (g in 1:nG)
di[g]=sum(abs(s1[g,]-s2[g,]))-abs(s1[g,g]-s2[g,g])
di=di/(nG-1)
perm.di=rep(0,nG)
count.ind=rep(0,nG)
for (i in 1:num.permutations){
i1=as.vector(sample(1:n,n1))
perm.X1=X[,i1]
perm.X2=X[,-i1]
perm.s1=plsnet(perm.X1)
perm.s2=plsnet(perm.X2)
for (g in 1:nG)
perm.di[g]=sum(abs(perm.s1[g,]-perm.s2[g,]))-
abs(perm.s1[g,g]-perm.s2[g,g])
perm.di=perm.di/(nG-1)
for (g in 1:nG)
if (perm.di[g]>=di[g])
count.ind[g]=count.ind[g]+1
}
p.value.ind=rep(0,nG)
for (g in 1:nG){
p.value.ind[g]=count.ind[g]/num.permutations
}
list(pp=p.value.ind,d=di)
}