#Comparing Modularity Models: 8 group power test
library(geomorph)
library(EMMLi)
source('compare.CR.r')

#simulate data with known covariance
sim.data<-function(n,s){
  a<-matrix(rnorm(n*nrow(s)),ncol=nrow(s))
  b<-a%*%chol(s)
  return(b)
}

#builds with Sw and Sb with random noise in all components
s.build<-function(p,sb,sw,mod.gp){
  s<-matrix(sb,nrow=p, ncol=p)
  for (i in 1:nlevels(mod.gp)){
    s[which(mod.gp==levels(mod.gp)[i]),which(mod.gp==levels(mod.gp)[i])]<-sw
  }
  low<-s[lower.tri(s)]
  for (i in 1:length(low)){
    low[i]<-rnorm(1,mean=low[i],sd=0.01)  
  }
  s[lower.tri(s)]<-low
  s[upper.tri(s)]<-t(s)[upper.tri(s)] 
  diag(s)<-rnorm(p,mean=1,sd=0.01)
  s<-matrix(Matrix::nearPD(s)$mat,nrow=p,ncol=p)
  return(s)
}


#FOR random covariance matrix: 8 module
iter<-999
N <- c(50, 100, 250)
p <- c(16,32,64)
nsimul <- 500
###### Power Simulation
s.b <- 0.6 
s.w <- c(0,.025,.05,.075,.1) #NOTE, this is ADDED to SB

CR.Power <- CR.Misspec <- EMMLi.Power <- EMMLi.Misspec <- array(NA, dim = c(length(p), length(N),length(s.w)))
row.names(CR.Power) <- row.names(CR.Misspec) <- row.names(EMMLi.Power) <- row.names(EMMLi.Misspec) <- paste("p_",p,sep="")
colnames(CR.Power) <- row.names(CR.Misspec) <-colnames(EMMLi.Power) <- colnames(EMMLi.Misspec) <- paste("N_",N,sep="")
dimnames(CR.Power)[[3]] <- dimnames(CR.Misspec)[[3]] <- dimnames(EMMLi.Power)[[3]] <- dimnames(EMMLi.Misspec)[[3]] <- paste("SW_",s.w,sep="")
EMMLi.best <- NULL
for (i in 1:length(p)){
  for (j in 1:length(N)){
    #hypotheses
    m2 <- gl(2,p[i]/2)
    m4 <- gl(4,p[i]/4)
    m8 <- gl(8,p[i]/8)
    
    #data 
    for (ii in 1:length(s.w)){
      S <- s.build(p[i],sb=s.b,sw=(s.b+s.w[ii]),mod.gp=m8)
      Y <- lapply(1:nsimul, function(k) sim.data(n=N[j],S))
      
      CR.m2 <- lapply(1:nsimul, function(k) modularity.test(Y[[k]], partition.gp = m2, iter = iter, print.progress = FALSE))
      CR.m4 <- lapply(1:nsimul, function(k) modularity.test(Y[[k]], partition.gp = m4, iter = iter, print.progress = FALSE))
      CR.m8 <- lapply(1:nsimul, function(k) modularity.test(Y[[k]], partition.gp = m8, iter = iter, print.progress = FALSE))
      CR.comp <- lapply(1:nsimul, function(k) compare.CR(CR.null = TRUE, CR.m2[[k]],CR.m4[[k]],CR.m8[[k]]))
      
      #Power: true (2-gp) vs. null
      CR.Power[i,j,ii] <- length(which(lapply(1:nsimul, function(k) CR.comp[[k]]$pairwise.P[1,4]<=0.05 )=="TRUE")) / nsimul  #how often is true model better than null?
      #Misspecification: best vs. true (8-gp)
      CR.min <- unlist(lapply(1:nsimul, function(k) which.min(CR.comp[[k]]$sample.z) ))
      CR.p <- unlist(lapply(1:nsimul, function(k) CR.comp[[k]]$pairwise.P[4,CR.min[k]]))
      CR.Misspec[i,j,ii] <- length(which(CR.p <=0.05)) / nsimul  #significant that are not input model
      
      ## EMMLi. best model
      X.1 <-as.character(seq(1:p[i]))
      models <- data.frame(X.1,m2, m4, m8)
      R <- lapply(1:nsimul, function(k) cor(Y[[k]]))
      EMMLi.out <- lapply(1:nsimul, function(k)  EMMLi(R[[k]],N[j], models)$results )
      best <- unlist(lapply(1:nsimul, function(k) which.min(EMMLi.out[[k]][,5])))
      #NOTE: for power, ANY model with 8 modules counted as best (a favorable view for EMMLi. I could be more strict)
      EMMLi.Power[i,j,ii] <- ifelse(ii==1,(1-(length(which(best== 1)) / nsimul)) ,length(which(best==8 | best==9 | best==10 | best==11)) / nsimul )
      EMMLi.Misspec[i,j,ii] <- 1 - EMMLi.Power[i,j,ii] #Misspecification is the converse
      
      best.out <- rep(0,11); best.out[sort(unique(best))] <- table(best)
      EMMLi.best <- rbind(EMMLi.best, c(p[i],N[j],s.w[ii], best.out))   #to ensure that all options are saved
      print(paste("i =",i))
      print(paste("j =",j))
      print(paste("ii =",ii))
    }
  }
}
colnames(EMMLi.best) <- c("p","N","s.w","1","2","3","4","5","6","7","8","9","10","11")

write.csv(EMMLi.best,"emmli.best-8gp-.6sb.csv")
write.csv(EMMLi.Power,"emmli.power-8gp-.6sb.csv")
write.csv(EMMLi.Misspec,"emmli.misspec-8gp-.6sb.csv")
write.csv(CR.Power, "CR.power-8gp-.6sb.csv")
write.csv(CR.Misspec,"CR.misspec-8gp-.6sb.csv")