Code for simulating figures 2 and 3

source("Functions For Simulations.R")

set.seed(12345)

n <- 500 #total number of vertices in one graph
kshuff <- c(10, 50, 100, 150, 200) #number shuffled in null
lshuff <- c(1, 5, 10, 25, 50, 100, 150, 200) #number shuffled in alt
prob1 <- matrix(.55, 250, 250) #block1 probabilities
prob2 <- matrix(.45, 250, 250) #block2 probabilities
prob3 <- matrix(.4, 250, 250) #block1 and 2 probabilities
P <- cbind(rbind(prob1, prob3), rbind(prob3, prob2)) #block probability matrix

block1 <- c(1:250) #list of vertices in block 1
block2 <- c(251:500) #list of vertices in block 2
nMC<-200 #total number of repetitions

null_reg <- list()
alt_reg <- list()
null_adj <- list()
alt_adj <- list()

error<-c(0.01, 0.02, 0.03, -0.01, -0.02, -0.03) #error levels
kk<-length(kshuff)
for(i in 1:length(error))
{
  null_reg[[i]]<-matrix(0,nMC,kk)
  null_adj[[i]]<-matrix(0,nMC,kk)
  
  alt_reg[[i]]<-list()
  alt_reg[[i]][[1]]<-matrix(0,nMC,3)
  alt_reg[[i]][[2]]<-matrix(0,nMC,5)
  alt_reg[[i]][[3]]<-matrix(0,nMC,6)
  alt_reg[[i]][[4]]<-matrix(0,nMC,7)
  alt_reg[[i]][[5]]<-matrix(0,nMC,8)
  
  alt_adj[[i]]<-list()
  alt_adj[[i]][[1]]<-matrix(0,nMC,3)
  alt_adj[[i]][[2]]<-matrix(0,nMC,5)
  alt_adj[[i]][[3]]<-matrix(0,nMC,6)
  alt_adj[[i]][[4]]<-matrix(0,nMC,7)
  alt_adj[[i]][[5]]<-matrix(0,nMC,8)
}

for (i in 1:length(error))
{
  for(h in 1:nMC)
  {
    for(j in 1:length(kshuff))
    {
      # make the graphs
      X <- sample_correlated_ieg_pair(n, P, P)$graph1
      A <- X[]
      Xase <- ase(A,2)$X
      P1hat<-Xase%*%t(Xase)
      P1hat[P1hat<0]<-0
      P1hat[P1hat>1]<-1
      
      X2 <- sample_correlated_ieg_pair(n, P, P)$graph1
      B <- X2[]
      X2ase <- ase(B,2)$X
      P2hat<-X2ase%*%t(X2ase)
      P2hat[P2hat<0]<-0
      P2hat[P2hat>1]<-1
      
      k<-kshuff[j]
      ord <- 1:n #original order of vertices
      b1_null <- sample(block1, k) #picking vert. in block 1 to shuff in H0
      b2_null <- sample(block2, k) #picking vert. in block 2 to shuff in H0
      ord_null <- shuff_ord(ord, b1_null, b2_null) #switching vertices for null
      P2hatshuff <- P2hat[ord_null,ord_null]
      
      null_reg[[i]][h,j] <- norm(as.matrix(P1hat - P2hatshuff), "F")^2
      null_adj[[i]][h,j] <- norm(as.matrix(A - B[ord_null, ord_null]), "F")^2
      
      s <- sum(lshuff<= k)
      
      for(q in 1:s)
      {
        E <- matrix(error[i], n, n)
        P_alt<- P + E
        X_alt <- sample_correlated_ieg_pair(n, P_alt,P_alt)$graph1
        C <- X_alt[]
        X_alt_ase <- ase(C,2)$X
        P3hat<-X_alt_ase%*%t(X_alt_ase)
        P3hat[P3hat<0]<-0
        P3hat[P3hat>1]<-1
        
        l<-lshuff[q]
        b1_alt <- sample(block1, l) 
        b2_alt <- sample(block2, l) 
        ord_alt <- shuff_ord(ord, b1_alt, b2_alt)
        P3hatshuff <- P3hat[ord_alt, ord_alt]
        
        alt_reg[[i]][[j]][h,q] <- norm(as.matrix(P1hat - P3hatshuff), "F")^2
        alt_adj[[i]][[j]][h,q] <- norm(as.matrix(A - C[ord_alt, ord_alt]), "F")^2
      }
    }
  }
}

#====================================================================================================================
#-----------------------------------------Calculating critical values for null---------------------------------------
#====================================================================================================================

null_crit_reg<-list()
null_crit_adj<-list()

for(i in 1:length(error))
{
  null_crit_reg[[i]]<-rep(0,5)
  null_crit_adj[[i]]<-rep(0,5)
  
  for(j in 1:length(kshuff))
  {
    null_crit_reg[[i]][j]<-quantile(null_reg[[i]][,j], 0.95)
    null_crit_adj[[i]][j]<-quantile(null_adj[[i]][,j], 0.95)
  }
}

#=====================================================================================================================
#---------------------------Calculating power by counting percentage of entries in alt -------------------------------
#-------------------------that are greater than their respective critical values in null------------------------------
#=====================================================================================================================
Pow_reg <-list()
Pow_adj <-list()

# Compute least favorable shuffling under H0
for(i in 1:length(error)){
  for(j in 1:length(kshuff)){
    null_crit_reg[[i]][j]<-max(null_crit_reg[[i]][1:j])
    null_crit_adj[[i]][j]<-max(null_crit_adj[[i]][1:j])
  }
}

for(i in 1:length(error))
{
  Pow_reg[[i]]<-matrix(NA,5,8)
  Pow_adj[[i]]<-matrix(NA,5,8)
  
  for(j in 1:length(kshuff))
  {
    kk<-kshuff[j]
    ll<-sum(lshuff<=kk)
    
    for(q in 1:ll)
    {
      Pow_reg[[i]][j,q]<-sum(alt_reg[[i]][[j]][,q] > null_crit_reg[[i]][j])/nMC
      Pow_adj[[i]][j,q]<-sum(alt_adj[[i]][[j]][,q] > null_crit_adj[[i]][j])/nMC
    }
  }
}
Pow_reg_sd<-list()
Pow_adj_sd<-list()

for(i in 1:length(error)){
  Pow_reg_sd[[i]] <- matrix(0,length(kshuff), length(lshuff))
  Pow_adj_sd[[i]] <- matrix(0,length(kshuff), length(lshuff))
  for(j in 1:length(kshuff))
  {
    kk<-kshuff[j]
    ll<-sum(lshuff<=kk)
    
    for(q in 1:ll)
    {
      m1<-sum(alt_reg[[i]][[j]][,q] > null_crit_reg[[i]][j])/nMC
      Pow_reg_sd[[i]][j,q]<-sqrt(m1*(1-m1)/nMC)
      m2<-sum(alt_adj[[i]][[j]][,q] > null_crit_adj[[i]][j])/nMC
      Pow_adj_sd[[i]][j,q]<-sqrt(m2*(1-m2)/nMC)
    }
  }
}
#===================================================================================================================
#-------------------------------------------Plotting Power-------------------------------------------------------
#===================================================================================================================

library(ggplot2)
library(reshape2)
library(viridis)
library(colorspace)
library(extrafont)
this <- rev(qualitative_hcl(5, "Dark 3"))
t <- "Error ="
t_adj <- "Error = " 

plot <- list()
plotadj <- list()
for (i in 1:length(Pow_reg))
{
  p <- Pow_reg[[i]]
  pp <- Pow_reg_sd[[i]]
  e1 <- as.data.frame(t(rbind(2*lshuff,p)))
  e2 <- as.data.frame(t(rbind(2*lshuff,pp)))
  colnames(e1) <- c("lshuff","20", "100", "200", "300", "400")
  colnames(e2) <- c("lshuff","20", "100", "200", "300", "400")
  e11 <- melt(e1, id.vars = "lshuff")
  e22 <- melt(e2, id.vars = "lshuff")
  colnames(e22) <- c("lshuff","variable", "se")
  e11$se<-e22$se
  e11$variable <- factor(e11$variable, levels = unique(e11$variable))
  e <- error[i]
  tfinal <- paste(t, as.character(e), sep=" ")
  
  plot[[i]] <- ggplot(e11, aes(x=lshuff, y=value, color=variable)) + 
    geom_line(aes(linetype=variable),size=1.05) +  
    geom_point(aes(shape=variable),size=4)+ ylim(0,1) +
    geom_errorbar(aes(ymin=value-2*se,ymax=value+2*se, width=20))+
    labs(family="serif", title=tfinal, x="Number of Vertices Shuffled in H1", y="Power")   + 
    theme_classic() + 
    theme(plot.title = element_text(family="serif", face = "bold", size = 17, hjust=0.5), 
          panel.background = element_rect(fill = "#FFFFFF"), 
          legend.title = element_text(family="serif", face="bold", color = "#000000", size = 12),
          legend.background = element_rect(fill = "#FFFFFF"), 
          axis.title = element_text(family = "serif", size = 15, colour = "#000000"),
          plot.background = element_rect(fill = "#FFFFFF"),
          axis.line = element_line(colour = "#000000", size = 1), 
          legend.text = element_text( family="serif", color = "#000000", size=15)) +  
          scale_color_manual(values=c("20"="#B675E0","100"="#00A6CA", "200"="#00AA5A", "300"="#AA9000", "400"="#E16A86"),breaks = 
          levels(e11$variable), name="Number shuffled in H0")+
          scale_shape_manual(values=c("20"= 18,"100"=16, "200"= 8, "300"= 17, "400"= 15),breaks = levels(e11$variable), name="Number shuffled in H0")+
    scale_linetype_manual(values=c("20"= 1,"100"=2, "200"= 3, "300"= 4,
                                         "400"=5),breaks = levels(e11$variable), 
                                name="Number shuffled in H0")
  
  p_adj <- Pow_adj[[i]]
  pp_adj <- Pow_adj_sd[[i]]
  e1_adj <- as.data.frame(t(rbind(2*lshuff,p_adj)))
  e2_adj <- as.data.frame(t(rbind(2*lshuff,pp_adj)))
  colnames(e1_adj) <- c("lshuff", "20", "100", "200", "300", "400")
  colnames(e2_adj) <- c("lshuff", "20", "100", "200", "300", "400")
  e11_adj <- melt(e1_adj, id.vars = "lshuff")
  e22_adj <- melt(e2_adj, id.vars = "lshuff")
  colnames(e22_adj) <- c("lshuff","variable", "se")
  e11_adj$se<-e22_adj$se
  e11_adj$variable <- factor(e11_adj$variable, levels = unique(e11_adj$variable))
  t_adjfinal <- paste0(t_adj, as.character(error[i]))
  
  plotadj[[i]] <- ggplot(e11_adj, aes(x=lshuff, y=value, color=variable)) + 
    geom_line(aes(linetype=variable),size=1.05) +  
    geom_point(aes(shape=variable),size=4)+ ylim(0,1) +
    geom_errorbar(aes(ymin=value-2*se,ymax=value+2*se), width=20)+
    labs(family="serif", title=t_adjfinal, x="Number of Vertices Shuffled in H1", y="Power")   + 
    theme_classic() + 
    theme(plot.title = element_text(family="serif", face = "bold", size = 17, hjust=0.5), 
          panel.background = element_rect(fill = "#FFFFFF"), 
          legend.title = element_text(family="serif", face="bold", color = "#000000", size = 12),
          legend.background = element_rect(fill = "#FFFFFF"), 
          axis.title = element_text(family = "serif", size = 15, colour = "#000000"),
          plot.background = element_rect(fill = "#FFFFFF"),
          axis.line = element_line(colour = "#000000", size = 1), 
          legend.text = element_text( family="serif", color = "#000000", size=15)) +  
          scale_color_manual(values=c("20"="#B675E0","100"="#00A6CA", "200"="#00AA5A", "300"="#AA9000", "400"="#E16A86"),breaks = 
          levels(e11_adj$variable), name="Number shuffled in H0")+
          scale_shape_manual(values=c("20"= 18,"100"=16, "200"= 8, "300"= 17, "400"= 15),breaks = levels(e11_adj$variable), name="Number shuffled in H0") +
    scale_linetype_manual(values=c("20"= 1,"100"=2, "200"= 3, "300"= 4,
                                         "400"=5),breaks = levels(e11_adj$variable), 
                                name="Number shuffled in H0")
}

library(patchwork)
## 
## Attaching package: 'patchwork'
## The following object is masked from 'package:MASS':
## 
##     area
plot[[1]] + plot[[2]] + plot[[3]] + plot[[4]] + plot[[5]] + plot[[6]] + plot_layout(ncol = 3, guides =
"collect")

plotadj[[1]] + plotadj[[2]] + plotadj[[3]] + plotadj[[4]] + plotadj[[5]] + plotadj[[6]] + plot_layout(ncol = 3, guides =
"collect")