#=========================================
# Function definitions
#=========================================

# ================================================================================
#   Cohesiveness/Conformity/Factorizability etc
# ================================================================================

# ===================================================
# Check if adj is a valid adjacency matrix:  square matrix, non-negative entries, symmetric and no missing entries.
# Parameters: 
#   adj - the input adjacency matrix
#   tol - the tolerence level to measure the difference from 0 (symmetric matrix: upper diagonal minus lower diagonal)
# Remarks:
#   1. This function is not supposed to be used directly. Instead, it should appear in function definitions.
#   2. We release the requirement that the diagonal elements be 1 or 0. Users should assign appropriate values
#      at the beginning of their function definitions.
# Usage:
# if(!is.adjmat(adj)) stop("The input matrix is not a valid adjacency matrix!")

if( exists("is.adjmat") ) rm(is.adjmat);
is.adjmat = function(adj, tol=10^(-15)){
	n=dim(adj)
	is.adj=1
	if (n[1] != n[2]){ message("The adjacency matrix is not a square matrix!"); is.adj=0;}
	if ( sum(is.na(adj))>0 ){ message("There are missing values in the adjacency matrix!"); is.adj=0;}
	if ( sum(adj<0)>0 ){ message("There are negative entries in the adjacency matrix!"); is.adj=0;}
	if ( max(abs(adj-t(adj))) > tol){ message("The adjacency matrix is not symmetric!"); is.adj=0;}
			#if ( max(abs(diag(adj)-1)) > tol){ message("The diagonal elements are not all one!"); is.adj=0;}
			#The last criteria is removed because of different definitions on diagonals with other papers.
			#Always let "diagonal=1" INSIDE the function calls when using functions for Factorizability paper.
	is.adj
}


# ===================================================
# NPC.direct=function(adj)
# Calculates the square root of Normalized Product Connectivity (NPC), by way of definition of NPC. ( \sqrt{t})
# Parameters: 
#   adj - the input adjacency matrix
#   tol - the tolerence level to measure the difference from 0 (zero off-diagonal elements)
# Output:
#   v1 - vector, the square root of NPC
# Remarks:
#   1. The function requires that the off-diagonal elements of the adjacency matrix are all non-zero.
#   2. If any of the off-diagonal elements is zero, use the function NPC.iterate().
#   3. If the adjacency matrix is 2 by 2, then a warning message is issued and vector of sqrt(adj[1,2]) is returned.
#   4. If the adjacency matrix is a ZERO matrix, then a warning message is issued and vector of 0 is returned.

if( exists("NPC.direct") ) rm(NPC.direct);
NPC.direct=function(adj){
	if(!is.adjmat(adj)) stop("The input matrix is not a valid adjacency matrix!")
	n=dim(adj)[1]
	if(n==2) {
		warning("The adjacecny matrix is only 2 by 2. NPC may not be unique!")
		return(rep(sqrt(adj[1,2]),2))
	}
	diag(adj)=0
	if(!sum(adj>0)){
	  warning("The adjacency matrix is a ZERO matrix!")
	  return(rep(0,n))
	}
	diag(adj)=1
	if(sum(adj==0)) stop("There is zero off--diagonal element! Please use the function NPC.iterate().")
	log10.prod.vec=function(vec){
		prod=0
		for(i in 1:length(vec) )
			prod=prod+log10(vec[i])
		prod
	}
	off.diag=as.vector(as.dist(adj))
	prod1=log10.prod.vec(off.diag)
	v1=rep(-666, n)
	for(i in 1:n){
		prod2=prod1-log10.prod.vec(adj[i,])
		v1[i]=10^(prod1/(n-1)-prod2/(n-2))
	}
	v1
}

# ===================================================
# NPC.iterate=function(adj, loop=10^(10), tol=10^(-10))
# Calculates the square root of Normalized Product Connectivity, by way of iteration algorithm. ( \sqrt{t})
# Parameters: 
#   adj - the input adjacency matrix
#   loop - the maximum number of iterations before stopping the algorithm
#   tol - the tolerence level to measure the difference from 0 (zero off-diagonal elements)
# Output:
#   v1 - vector, the square root of NPC
#   loop - integer, the number of iterations taken before convergence criterion is met
#   diff - scaler, the maximum difference between the estimates of 'v1' in the last two iterations
# Remarks:
#   1. Whenever possible, use NPC.direct().
#   2. If the adjacency matrix is 2 by 2, then a warning message is issued.
#   3. If the adjacency matrix is a ZERO matrix, then a warning message is issued and vector of 0 is returned.

if( exists("NPC.iterate") ) rm(NPC.iterate);
NPC.iterate=function(adj, loop=10^(10), tol=10^(-10)){
	if(!is.adjmat(adj)) stop("The input matrix is not a valid adjacency matrix!")
	n=dim(adj)[1]
	if(n==2) warning("The adjacecny matrix is only 2 by 2. NPC may not be unique!")
	diag(adj)=0
	if(max(abs(adj))<tol){
	  warning("The adjacency matrix is a ZERO matrix!")
	  return(rep(0,n))
	}
	diff=1
	k=apply(adj, 2, sum)
	v1=k/sqrt(sum(k)) # first-step estimator for v-vector (initial value)
	i=0
	while( loop>i && diff>tol ){
	i=i+1
	diag(adj)=v1^2
	svd1=svd(adj) # Spectral Decomposition
	v2=sqrt(svd1$d[1])*abs(svd1$u[,1])
	diff=max(abs(v1-v2))
	v1=v2
	}
	list(v1=v1,loop=i,diff=diff)
}

# ===================================================
# FDADJ = function(adj)
# Calculate module specific quantities.
# The name is the abbreviation for "Factorizability Decomposition of ADJacency matrix"
# Parameters: 
#   adj - The input [adjacency] matrix; The general matrix M in the definition of network functions.
#         The diagonals of M may not be 0 and M may not be symmetric.
# Output: the output is a 'list' object containing the following objects:
#   Size - integer, the Module Size
#   Factorizability - scaler, the Factorizability of the module
#   eta1 - scaler, the Squared Mean Conformity
#   ClusterCoef - vector, the Clustering Coefficients
#   MConformity - vector, the Module Conformity
#   kVar - vector, the Within Module Connectivities Variance
#   kWithin - vector, the Within Module Connectivities
#   kEst - vector, (kApprox - NPC)
#   kApprox - vector, the approximated Within Module Connectivities by scaled conformity
#   summary - vector, summarized information of the module
#         "Size", "Factorizability", "Uniformity", "Uniformity.CF", "Tightness",
#	        "Tightness.CF", "C.CF", "ClusterCoef.Mean", "kWithin.Mean", "Corr(k, kCF)"
# Remarks:
#   1. In default, the function NPC.direct() is used to calculate NPC, with the definition of NPC.
#   2. When there are zero elements in the adjacency matrix, NPC.iterate() will be used to calculate NPC, an iterative algorithm. 
if( exists("FDADJ") ) rm(FDADJ);
FDADJ = function(adj) {
	if(!is.adjmat(adj)) stop("The input matrix is not a valid adjacency matrix!")
	diag(adj)=0 # Therefore adj=A-I.
	
	### Fundamental Network Concepts
	Size=dim(adj)[1]
	Connectivity=apply(adj, 2, sum) # Within Module Connectivities
	Density=sum(Connectivity)/(Size*(Size-1))
	Centralization=Size*(max(Connectivity)-mean(Connectivity))/((Size-1)*(Size-2))
	Heterogeneity=sqrt(Size*sum(Connectivity^2)/sum(Connectivity)^2-1)
	ClusterCoef=ClusterCoef.fun(adj)
	
	fMAR=function(v) sum(v^2)/sum(v)
	MAR=apply(adj, 1, fMAR)
	#CONNECTIVITY=Connectivity/max(Connectivity)
	
	### Conformity-Based Network Concepts	
	Conformity=NPC.iterate(adj)$v1
	Factorizability=1- sum( (adj-outer(Conformity,Conformity)+ diag(Conformity^2))^2 )/sum(adj^2)
	Connectivity.CF=sum(Conformity)*Conformity-Conformity^2
	Density.CF=sum(Connectivity.CF)/(Size*(Size-1))
	Centralization.CF=Size*(max(Connectivity.CF)-mean(Connectivity.CF))/((Size-1)*(Size-2))
	Heterogeneity.CF=sqrt(Size*sum(Connectivity.CF^2)/sum(Connectivity.CF)^2-1)
	#ClusterCoef.CF=ClusterCoef.fun(outer(Conformity,Conformity)-diag(Conformity^2) )
	ClusterCoef.CF=c(NA, Size)
  for(i in 1:Size )
    ClusterCoef.CF[i]=( sum(Conformity[-i]^2)^2 - sum(Conformity[-i]^4) )/( sum(Conformity[-i])^2 - sum(Conformity[-i]^2) )


	### Approximate Conformity-Based Network Concepts	
	Connectivity.CF.App=sum(Conformity)*Conformity
	Density.CF.App=sum(Connectivity.CF.App)/(Size*(Size-1))
	Centralization.CF.App=Size*(max(Connectivity.CF.App)-mean(Connectivity.CF.App))/((Size-1)*(Size-2))
	Heterogeneity.CF.App=sqrt(Size*sum(Connectivity.CF.App^2)/sum(Connectivity.CF.App)^2-1)
	ClusterCoef.CF.App=(sum(Conformity^2)/sum(Conformity))^2
	
	output=list(Conformity=Conformity, ClusterCoef=ClusterCoef, ClusterCoef.CF=ClusterCoef.CF,
	  kWithin=Connectivity, kWithin.CF=Connectivity.CF, kWithin.CF.App=Connectivity.CF.App, MAR=MAR)
	
	module=rep(NA,16)
	module.names=c("Size", "Factorizability",
	  "Density", "Density.CF", "Density.CF.App",
	  "Centralization", "Centralization.CF", "Centralization.CF.App",
	  "Heterogeneity", "Heterogeneity.CF", "Heterogeneity.CF.App",
		"ClusterCoef.Mean", "ClusterCoef.CF.Mean", "ClusterCoef.CF.App", "Conformity.Mean", "kWithin.Mean")
	module[1]=Size
	module[2]=Factorizability
	module[3]=Density
	module[4]=Density.CF
	module[5]=Density.CF.App
	module[6]=Centralization
	module[7]=Centralization.CF
	module[8]=Centralization.CF.App
	module[9]=Heterogeneity
	module[10]=Heterogeneity.CF
	module[11]=Heterogeneity.CF.App
	module[12]=mean(ClusterCoef)
	module[13]=mean(ClusterCoef.CF)
	module[14]=ClusterCoef.CF.App
	module[15]=mean(Conformity)
	module[16]=mean(Connectivity)

	output$summary=module
	output$summary.names=module.names
	output
}


# ================================================================================
# Additional code from Jun Dong:
#   Modified functions in "pairs()" function.
# ================================================================================


# ===================================================
# This function is used in "pairs()" function. 
# The problem of the original panel.cor is that when the correlation 
# coefficient is very small, the lower panel will have a large font 
# instead of a mini-font in a saved .ps file. This new function uses 
# a format for corr=0.2 when corr<0.2, but it still reports the original 
# value of corr, with a minimum format.

if(exists("panel.cor1")) rm(panel.cor1);
panel.cor1=function(x, y, digits=2, prefix="", cex.cor){
    usr <- par("usr"); on.exit(par(usr))
    par(usr = c(0, 1, 0, 1))
    tmp=!is.na(x) & !is.na(y) & x!=Inf & x!=-Inf & y!=Inf & y!=-Inf
    r <- abs(cor(x[tmp], y[tmp]))
    txt <- format(c(r, 0.123456789), digits=digits)[1]
    txt <- paste(prefix, txt, sep="")
    txt1=txt
    r1=r
    if (r<0.2) {
        r1=0.2
        txt1 <- format(c(r1, 0.123456789), digits=digits)[1]
        txt1 <- paste(prefix, txt1, sep="")
        }
    if(missing(cex.cor)) cex <- 0.8/strwidth(txt1)
    cex = cex * r1
    r <- round(r, digits)
    txt <- format(c(r, 0.123456789), digits=digits)[1]
    txt <- paste(prefix, txt, sep="")
    text(0.5, 0.5, txt, cex=cex)
}


# ===================================================
# This function is used in "pairs()" function.
# The only difference with the original panel.smooth is a
# smaller value of cex = 0.7 instead of 1.
if(exists("panel.smooth1")) rm(panel.smooth1);
panel.smooth1 = function (x, y, col = par("col"), bg = NA, pch = par("pch"), 
    cex = 0.7, col.smooth = "red", span = 2/3, iter = 3, ...) 
{
    points(x, y, pch = pch, col = col, bg = bg, cex = cex)
    ok <- is.finite(x) & is.finite(y)
    if (any(ok)) 
        lines(stats::lowess(x[ok], y[ok], f = span, iter = iter), 
            col = col.smooth, ...)
}
# ===================================================
# This function is used in "pairs()" function.
# It allows to fit a linear line in the scatter plot.

if(exists("panel.lm")) rm(panel.lm);
panel.lm = function (x, y, col = par("col"), bg = NA, pch = par("pch"), 
    cex = 0.7, col.lm = "red", ...) 
{
    points(x, y, pch = pch, col = col, bg = bg, cex = cex)
    ok <- is.finite(x) & is.finite(y)
    if (any(ok)) 
        abline(lm(y[ok] ~ x[ok])$coefficients, col = col.lm, ...)
}

# ================================================================================
#  Steve Horvath's code for Module Principal Component Analysis
# ================================================================================

# ===================================================
#The function ModulePrinComps2 finds the first principal component in every 
# module defined by the colors of the input vector "couleur"
# This requires the R library impute.
# New version of Email, 1-18-2005
if(exists("ModulePrinComps2")) rm(ModulePrinComps2);
ModulePrinComps2=function(datexpr,couleur) {
modlevels=levels(factor(couleur))
PrinComps=data.frame(matrix(666,nrow=dim(datexpr)[[1]],ncol= length(modlevels))) 
varexplained= data.frame(matrix(666,nrow= 5,ncol= length(modlevels)))
names(PrinComps)=paste("PC",modlevels,sep="")
for(i in c(1:length(modlevels)) ){
#print(i)
modulename    = modlevels[i]
 restrict1= as.character(couleur)== modulename
 # in the following, rows are genes and columns are samples
    datModule=t(datexpr[, restrict1])
    datModule=impute.knn(as.matrix(datModule))
datModule=t(scale(t(datModule)))
svd1=svd(datModule)
    mtitle=paste("PCs of ", modulename," module", sep="")
varexplained[,i]= (svd1$d[1:5])^2/sum(svd1$d^2)
# this is the first principal component
    pc1=svd1$v[,1]
 signh1=sign(sum(cor(pc1,  t(datModule))))
 if (signh1 != 0)  pc1=signh1* pc1
PrinComps[,i]= pc1
}
list(PrinComps=PrinComps, varexplained=varexplained)
}


# ===========================================================
# The function PickSoftThreshold allows one to estimate the power parameter when using
# a soft thresholding approach with the use of the power function AF(s)=s^Power
# The function PickSoftThreshold allows one to estimate the power parameter when using
# a soft thresholding approach with the use of the power function AF(s)=s^Power
# The removeFirst option removes the first point (k=1, P(k=1)) from the regression fit.
if (exists("PickSoftThreshold1")) rm(PickSoftThreshold1);
PickSoftThreshold1=function(cormat,RsquaredCut=0.85, powervector=c(seq(1,10,by=1),seq(12,20,by=2)) ) {
removeFirst=FALSE
cormat=abs(cormat)
no.genes <- dim(cormat)[[1]]
colname1=c("Power", "scale law R^2" ,"slope", "truncated R^2","mean(k)","median(k)","max(k)")
datout=data.frame(matrix(666,nrow=length(powervector),ncol=length(colname1) ))
names(datout)=colname1
datout[,1]=powervector
if(exists("fun1")) rm(fun1)
fun1=function(corx) {
out1=rep(NA, length(powervector) )
for (j in c(1:length(powervector))) {out1[j]=sum(corx^powervector[j])}
out1
} # end of fun1
datk=t(apply(cormat,2,fun1))
for (i in c(1:length(powervector) ) ){
nolinkshelp <- datk[,i]-1
no.breaks=10
cut2=cut(nolinkshelp,no.breaks)
binned.k=tapply(nolinkshelp,cut2,mean)
freq1=as.vector(tapply(nolinkshelp,cut2,length)/length(nolinkshelp))
xx= as.vector(log10(binned.k))
if(removeFirst) {freq1=freq1[-1]; xx=xx[-1]}
plot(xx,log10(freq1+.000000001),xlab="log10(k)",ylab="log10(p(k))" )
lm1= lm(as.numeric(log10(freq1+.000000001))~ xx )
lm2=lm(as.numeric(log10(freq1+.000000001))~ xx+I(10^xx) )
datout[i,2]=summary(lm1)$adj.r.squared 
datout[i,3]=summary(lm1)$coefficients[2,1]  
datout[i,4]=summary(lm2)$adj.r.squared
datout[i,5]=mean(nolinkshelp)
datout[i,6]= median(nolinkshelp)
datout[i,7]= max(nolinkshelp) 
} 
datout=signif(datout,3) 
print(data.frame(datout));
# the cut-off is chosen as smallest cut with R^2>RsquaredCut 
ind1=datout[,2]>RsquaredCut
indcut=NA
indcut=ifelse(sum(ind1)>0,min(c(1:length(ind1))[ind1]),indcut)
# this estimates the power value that should be used. 
# Don't trust it. You need to consider slope and mean connectivity as well!
power.estimate=powervector[indcut][[1]]
list(power.estimate, data.frame(datout));
}

# ===================================================
#The function TOMdist1 computes a dissimilarity 
# based on the topological overlap matrix (Ravasz et al)
# Input: an Adjacency matrix with entries in [0,1]
# This function is different with TOMdist1 by not assigning the diagonals as 1.
if(exists("TOMdist2")) rm(TOMdist2);
TOMdist2=function(adjmat1, maxADJ=FALSE) {
diag(adjmat1)=0;
adjmat1[is.na(adjmat1)]=0;
maxh1=max(as.dist(adjmat1) ); minh1=min(as.dist(adjmat1) ); 
if (maxh1>1 | minh1 < 0 ) {print(paste("ERROR: the adjacency matrix contains entries that are larger than 1 or smaller than 0!!!, max=",maxh1,", min=",minh1)) } else { 
if (  max(c(as.dist(abs(adjmat1-t(adjmat1)))))>0   ) {print("ERROR: non-symmetric adjacency matrix!!!") } else { 
kk=apply(adjmat1,2,sum)
maxADJconst=1
if (maxADJ==TRUE) maxADJconst=max(c(as.dist(adjmat1 ))) 
Dhelp1=matrix(kk,ncol=length(kk),nrow=length(kk))
denomTOM= pmin((Dhelp1),(t(Dhelp1)))   +(maxADJconst-adjmat1); 
gc();gc();
numTOM=(adjmat1 %*% adjmat1 +adjmat1);
#TOMmatrix=numTOM/denomTOM
# this turns the TOM matrix into a dissimilarity 
out1=1-as.matrix(numTOM/denomTOM) 
#diag(out1)=1
out1
}}
}

# ===================================================
# The function ModuleEnrichment2 creates a bar plot that shows whether modules are enriched with
# significant genes.
# More specifically, it reports the mean gene significance for each module.
# The gene significance can be a binary variable or a quantitative variable. 
# It also plots the 95% confidence interval of the mean (CI=mean +/- 1.96* standard error).
# It also reports a Kruskal Wallis P-value.
if( exists("ModuleEnrichment2") ) rm(ModuleEnrichment2);
ModuleEnrichment2=function(genesignif1,couleur,title1="gene significance across modules",labely="Gene Significance",boxplot=F,
cex.lab1=1, cex.main1=1) {
if (length(genesignif1) != length(couleur) ) print("Error: vectors don¡¯t have the same lengths") else {
if (boxplot != TRUE) {
mean1=function(x) mean(x,na.rm=T) 
means1=as.vector(tapply(genesignif1,couleur,mean1));
se1= as.vector(tapply(genesignif1,couleur,stderr1))
par(mfrow=c(1,1))
barplot(means1, names.arg=names(table(couleur) ),col= names(table(couleur) ), cex.lab=cex.lab1)
err.bp(as.vector(means1), as.vector(1.96*se1), two.side=T)} else {
boxplot(split(genesignif1,couleur),notch=T,varwidth=T, col= names(table(couleur) ))}

title(paste(title1,", p-value=", signif(kruskal.test(genesignif1,factor(couleur))$p.value,2)), cex.main=cex.main1,ylab=labely, cex.lab=cex.lab1)
}
} # end of function


# ===================================================
# NetworkConcepts = function(datE, power=1, Trait=NULL, SignedNetwork=F){
# datE: expression profiles with rows=samples and cols=genes/probesets
# power: for contruction of the weighted network
# Trait: the quantitative external trait
# SignedNetwork: an indicator of whether corr= 0.5 is different with corr=-0.5
if( exists("NetworkConcepts") ) rm(NetworkConcepts);
NetworkConcepts = function(datE, power=1, Trait=NULL, SignedNetwork=F) {
  if(SignedNetwork){
  	adj <- abs((cor(datE,use="p")+1)/2)^power
  } else{
    adj <- abs(cor(datE,use="p"))^power
  }
	diag(adj)=0 # Therefore adj=A-I.
	
	### Fundamental Network Concepts
	Size=dim(adj)[1]
	Connectivity=apply(adj, 2, sum) # Within Module Connectivities
	Density=sum(Connectivity)/(Size*(Size-1))
	Centralization=Size*(max(Connectivity)-mean(Connectivity))/((Size-1)*(Size-2))
	Heterogeneity=sqrt(Size*sum(Connectivity^2)/sum(Connectivity)^2-1)
	ClusterCoef=ClusterCoef.fun(adj)
	
	fMAR=function(v) sum(v^2)/sum(v)
	MAR=apply(adj, 1, fMAR)
	#CONNECTIVITY=Connectivity/max(Connectivity)
	
	### Conformity-Based Network Concepts	
	### Dong J, Horvath S (2007) Understanding Network Concepts in Modules, BMC Systems Biology 2007, 1:24
	Conformity=NPC.iterate(adj)$v1
	Factorizability=1- sum( (adj-outer(Conformity,Conformity)+ diag(Conformity^2))^2 )/sum(adj^2)
	Connectivity.CF=sum(Conformity)*Conformity-Conformity^2
	Density.CF=sum(Connectivity.CF)/(Size*(Size-1))
	Centralization.CF=Size*(max(Connectivity.CF)-mean(Connectivity.CF))/((Size-1)*(Size-2))
	Heterogeneity.CF=sqrt(Size*sum(Connectivity.CF^2)/sum(Connectivity.CF)^2-1)
	#ClusterCoef.CF=ClusterCoef.fun(outer(Conformity,Conformity)-diag(Conformity^2) )
	ClusterCoef.CF=c(NA, Size)
  for(i in 1:Size )
    ClusterCoef.CF[i]=( sum(Conformity[-i]^2)^2 - sum(Conformity[-i]^4) )/( sum(Conformity[-i])^2 - sum(Conformity[-i]^2) )

	### Approximate Conformity-Based Network Concepts	
	Connectivity.CF.App=sum(Conformity)*Conformity
	Density.CF.App=sum(Connectivity.CF.App)/(Size*(Size-1))
	Centralization.CF.App=Size*(max(Connectivity.CF.App)-mean(Connectivity.CF.App))/((Size-1)*(Size-2))
	Heterogeneity.CF.App=sqrt(Size*sum(Connectivity.CF.App^2)/sum(Connectivity.CF.App)^2-1)
	ClusterCoef.CF.App=(sum(Conformity^2)/sum(Conformity))^2

  ### Eigengene-based Network Concepts
  m1=ModulePrinComps1(datE)
  ConformityE=cor(datE,m1[[1]][,1],use="pairwise.complete.obs"); ConformityE=abs(ConformityE)^power; # Weighted Expression Conformity
  ConnectivityE=sum(ConformityE)*ConformityE; #Expression Connectivity
  DensityE=sum(ConnectivityE)/(Size*(Size-1)); #Expression Density
  CentralizationE=Size*(max(ConnectivityE)-mean(ConnectivityE))/((Size-1)*(Size-2)); #Expression Centralization
  HeterogeneityE=sqrt(Size*sum(ConnectivityE^2)/sum(ConnectivityE)^2-1); #Expression Heterogeneity
  ClusterCoefE=(sum(ConformityE^2)/sum(ConformityE))^2; ##Expression ClusterCoef
  MARE=ConformityE* sum(ConformityE^2)/sum(ConformityE)

  ### Significance measure only when Trait is available.
  if(!is.null(Trait)){
    EigengeneSignificance= abs(cor(Trait, m1[[1]], use="pairwise.complete.obs") )^power; EigengeneSignificance=EigengeneSignificance[1]
    GS= abs(cor(datE, Trait, use="pairwise.complete.obs") )^power; GS=GS[,1]
    GSE=ConformityE * EigengeneSignificance; GSE=GSE[,1]
    ModuleSignificance=mean(GS)
    ModuleSignificanceE=mean(GSE)
    K=Connectivity/max(Connectivity)
    HubGeneSignificance=sum(GS*K)/sum(K^2)
    KE=ConnectivityE/max(ConnectivityE)
    HubGeneSignificanceE= sum(GSE*KE)/sum(KE^2)
  }

  Summary=cbind(
    c(Density, Centralization, Heterogeneity, mean(ClusterCoef), mean(Connectivity)),
    c(DensityE, CentralizationE, HeterogeneityE, mean(ClusterCoefE), mean(ConnectivityE)),
    c(Density.CF, Centralization.CF, Heterogeneity.CF, mean(ClusterCoef.CF), mean(Connectivity.CF)),
    c(Density.CF.App, Centralization.CF.App, Heterogeneity.CF.App, mean(ClusterCoef.CF.App), mean(Connectivity.CF.App) ) )
  colnames(Summary)=c("Fundamental", "Eigengene-based", "Conformity-Based", "Approximate Conformity-based")
  rownames(Summary)=c("Density", "Centralization", "Heterogeneity", "Mean ClusterCoef", "Mean Connectivity")

  output=list(Summary=Summary, Size=Size, Factorizability=Factorizability, Eigengene=m1[[1]], VarExplained=m1[[2]][,1], Conformity=Conformity, ClusterCoef=ClusterCoef, Connectivity=Connectivity, MAR=MAR, ConformityE=ConformityE)
  if(!is.null(Trait)){
    output$GS=GS; output$GSE=GSE;
    Significance=cbind(c(ModuleSignificance, HubGeneSignificance, EigengeneSignificance), c(ModuleSignificanceE, HubGeneSignificanceE, NA))
    colnames(Significance)=c("Fundamental", "Eigengene-based")
    rownames(Significance)=c("ModuleSignificance", "HubGeneSignificance", "EigengeneSignificance")
    output$Significance=Significance
  }
	output
}