# Homework 2 for the University of Tulsa's CS-7863 Network Theory Course
# Subgraph Centrality Comparisons, Microstate Computations, and Entropy
# Professor: Dr. McKinney, Spring 2022
# Noah Schrick - 1492657

# Imports
#install.packages("igraph")
#install.packages("igraphdata")
#install.packages("reshape2")
library(igraph)
library(igraphdata)
library(reshape2)
data(karate)
data(kite)
source("katz_centrality.R")
source("self_estrada.R")

# 3 Networks: Karate, Kite, and Fig. 1a of the subgraph centrality paper
g.one <- karate
g.one.netname <- "Karate"

g.two <- kite
g.two.netname <- "Kite"

g.fig1a <- graph.ring(8)
for (idx in V(g.fig1a)) {V(g.fig1a)[idx]$name <- idx}
g.fig1a <- g.fig1a %>% add_edges(c(2,8, 3,6, 4,7, 1,5))
g.fig1a.netname <- "Fig1a"

####################### Part 1: Centrality Comparisons #########################
# Container to hold results for each centrality measure
centralities <- matrix(list(), nrow=3, ncol=5)
rownames(centralities) <- c(g.one.netname, g.two.netname, g.fig1a.netname)
colnames(centralities) <- c("Eigenvector", "Subgraph", "Betweenness", "Katz", 
                            "Subgraph-like Katz")

# Eigenvector - using igraph
centralities[[1,1]] <- eigen_centrality(g.one)$vector %>% sort(decreasing=T)
centralities[[2,1]] <- eigen_centrality(g.two)$vector %>% sort(decreasing=T)
centralities[[3,1]] <- eigen_centrality(g.fig1a)$vector %>% sort(decreasing=T)

# Subgraph - using igraph
centralities[[1,2]] <- subgraph.centrality(g.one) %>% sort(decreasing=T)
centralities[[2,2]] <- subgraph.centrality(g.two) %>% sort(decreasing=T)
centralities[[3,2]] <- subgraph.centrality(g.fig1a) %>% sort(decreasing=T)

# Betweenness - using igraph
centralities[[1,3]] <- betweenness(g.one) %>% sort(decreasing=T)
centralities[[2,3]] <- betweenness(g.two) %>% sort(decreasing=T)
centralities[[3,3]] <- betweenness(g.fig1a) %>% sort(decreasing=T)

# Katz
centralities[[1,4]] <- katz.cent(g.one) %>% sort(decreasing=T)
centralities[[2,4]] <- katz.cent(g.two) %>% sort(decreasing=T)
centralities[[3,4]] <- katz.cent(g.fig1a) %>% sort(decreasing=T)


# "Subgraph-like Katz"
centralities[[1,5]] <- sg.katz(g.one) %>% sort(decreasing=T)
centralities[[2,5]] <- sg.katz(g.two) %>% sort(decreasing=T)
centralities[[3,5]] <- sg.katz(g.fig1a) %>% sort(decreasing=T)


####################### Part 2: Microstates and Entropy ########################
# Create 3 100-node Erdos-Renyi random graph: 10%, 50%, and 100% attachment
# probabilities. Will use a Beta of 0.01, 0.5, and 1.0 for all 3 networks.
g.r.one <- erdos.renyi.game(100, .1)
# for (idx in V(g.r.one)) {V(g.r.one)[idx]$name <- idx}
g.r.one.netname <- "10% Attachment Probability"

g.r.two <- erdos.renyi.game(100, .5)
for (idx in V(g.r.two)) {V(g.r.two)[idx]$name <- idx}
g.r.two.netname <- "50% Attachment Probability"

g.r.three <- erdos.renyi.game(100, 1)
for (idx in V(g.r.three)) {V(g.r.three)[idx]$name <- idx}
g.r.three.netname <- "100% Attachment Probability"

# Container to hold results for each network
res <- matrix(list(), nrow=3, ncol=4)
rownames(res) <- c(g.r.one.netname, g.r.two.netname, g.r.three.netname)
colnames(res) <- c("Degree Distribution", "Estrada Index",
                            "Microstates p_i's", "Entropy")

# Degree Distribution-Not Dependent on Beta, but treating it for easier results
res[[1,1]] <- list("0.01"=degree_distribution(g.r.one) %>% sort(decreasing=T), 
                   "0.5"=degree_distribution(g.r.one) %>% sort(decreasing=T),
                   "1.0"=degree_distribution(g.r.one) %>% sort(decreasing=T))
res[[2,1]] <- list("0.01"=degree_distribution(g.r.two) %>% sort(decreasing=T), 
                   "0.5"=degree_distribution(g.r.two) %>% sort(decreasing=T),
                   "1.0"=degree_distribution(g.r.two) %>% sort(decreasing=T))
res[[3,1]] <- list("0.01"=degree_distribution(g.r.three) %>% sort(decreasing=T),
                   "0.5"=degree_distribution(g.r.three) %>% sort(decreasing=T), 
                   "1.0"=degree_distribution(g.r.three) %>% sort(decreasing=T))

# Compute EE
res[[1,2]] <- list("0.01"=estrada.index(g.r.one, 0.01), 
                   "0.5"=estrada.index(g.r.one, 0.5),
                   "1.0"=estrada.index(g.r.one, 1))
res[[2,2]] <- list("0.01"=estrada.index(g.r.two, 0.01),
                   "0.5"=estrada.index(g.r.two, 0.5),
                   "1.0"=estrada.index(g.r.two, 1))
res[[3,2]] <- list("0.01"=estrada.index(g.r.three, 0.01),
                   "0.5"=estrada.index(g.r.three, 0.5),
                   "1.0"=estrada.index(g.r.three, 1))

# Compute Microstates
res[[1,3]] <- list("0.01"=microstate.prob(g.r.one, 0.01), 
                   "0.5"=microstate.prob(g.r.one, 0.5),
                   "1.0"=microstate.prob(g.r.one, 1))
res[[2,3]] <- list("0.01"=microstate.prob(g.r.two, 0.01),
                   "0.5"=microstate.prob(g.r.two, 0.5),
                   "1.0"=microstate.prob(g.r.two, 1))
res[[3,3]] <- list("0.01"=microstate.prob(g.r.three, 0.01),
                   "0.5"=microstate.prob(g.r.three, 0.5),
                   "1.0"=microstate.prob(g.r.three, 1))

# Histogram Display
# 10% Attachment Probability Network
par(mfrow=c(3,1))
g.r<- unlist(res[[1,3]]["0.01"], use.names=FALSE)
hist(g.r,
    main="Microstate Histogram for a a 100-node Erdos-Renyi random
    graph with attachment probability 10% and Beta=0.01",
    xlab="Occupation Probability")
g.r <- unlist(res[[1,3]]["0.5"], use.names=FALSE)
hist(g.r,
     main="Microstate Histogram for a a 100-node Erdos-Renyi random
    graph with attachment probability 10% and Beta=0.5",
     xlab="Occupation Probability")
g.r <- unlist(res[[1,3]]["1.0"], use.names=FALSE)
hist(g.r,
     main="Microstate Histogram for a a 100-node Erdos-Renyi random
    graph with attachment probability 10% and Beta=1.0",
     xlab="Occupation Probability")
# 50% Attachment Probability Network
par(mfrow=c(3,1))
g.r<- unlist(res[[2,3]]["0.01"], use.names=FALSE)
hist(g.r,
     main="Microstate Histogram for a a 100-node Erdos-Renyi random
    graph with attachment probability 50% and Beta=0.01",
     xlab="Occupation Probability")
g.r <- unlist(res[[2,3]]["0.5"], use.names=FALSE)
hist(g.r,
     main="Microstate Histogram for a a 100-node Erdos-Renyi random
    graph with attachment probability 50% and Beta=0.5",
     xlab="Occupation Probability")
g.r <- unlist(res[[2,3]]["1.0"], use.names=FALSE)
hist(g.r,
     main="Microstate Histogram for a a 100-node Erdos-Renyi random
    graph with attachment probability 50% and Beta=1.0",
     xlab="Occupation Probability")
# 100% Attachment Probability Network
par(mfrow=c(3,1))
g.r<- unlist(res[[3,3]]["0.01"], use.names=FALSE)
hist(g.r,
     main="Microstate Histogram for a 100-node Erdos-Renyi random
    graph with attachment probability 100% and Beta=0.01",
     xlab="Occupation Probability")
g.r <- unlist(res[[3,3]]["0.5"], use.names=FALSE)
hist(g.r,
     main="Microstate Histogram for a 100-node Erdos-Renyi random
    graph with attachment probability 100% and Beta=0.5",
     xlab="Occupation Probability")
g.r <- unlist(res[[3,3]]["1.0"], use.names=FALSE)
hist(g.r,
     main="Microstate Histogram for a 100-node Erdos-Renyi random
    graph with attachment probability 100% and Beta=1.0",
     xlab="Occupation Probability")

# Entropy
res[[1,4]] <- list("0.01"=entropy(g.r.one, 0.01), 
                   "0.5"=entropy(g.r.one, 0.5),
                   "1.0"=entropy(g.r.one, 1))
res[[2,4]] <- list("0.01"=entropy(g.r.two, 0.01),
                   "0.5"=entropy(g.r.two, 0.5),
                   "1.0"=entropy(g.r.two, 1))
res[[3,4]] <- list("0.01"=entropy(g.r.three, 0.01),
                   "0.5"=entropy(g.r.three, 0.5),
                   "1.0"=entropy(g.r.three, 1))
y1 <- c(unlist(res[[1,4]]["0.01"], use.names=FALSE),
                unlist(res[[2,4]]["0.01"], use.names=FALSE),
                unlist(res[[3,4]]["0.01"], use.names=FALSE))

y2 <- c(unlist(res[[1,4]]["0.5"], use.names=FALSE),
                unlist(res[[2,4]]["0.5"], use.names=FALSE),
                unlist(res[[3,4]]["0.5"], use.names=FALSE))

y3 <-c(unlist(res[[1,4]]["1.0"], use.names=FALSE),
                 unlist(res[[2,4]]["1.0"], use.names=FALSE),
                 unlist(res[[3,4]]["1.0"], use.names=FALSE))

par(mfrow=c(3,1))

plot(x=c(0.1,0.5,1.0), y=y1, col="red", type="o", pch="o",
        xlab="Attachment Probability", ylab="Entropy",
        main="Entropies for a 100-node Erdos-Renyi random
    graph with varying attachment probability and Beta=0.01", lty=1,
     ylim=c(min(y1),max(y1)))

plot(x=c(0.1,0.5,1.0), y=y2, col="red", type="o", pch="o",
     xlab="Attachment Probability", ylab="Entropy",
     main="Entropies for a 100-node Erdos-Renyi random
    graph with varying attachment probability and Beta=0.5", lty=1, ylim=c(min(y2),max(y2)))

plot(x=c(0.1,0.5,1.0), y=y3, col="red", type="o", pch="o",
     xlab="Attachment Probability", ylab="Entropy",
     main="Entropies for a 100-node Erdos-Renyi random
    graph with varying attachment probability and Beta=1.0", lty=1, ylim=c(min(y3),max(y3)))


# igraph network entropy
g1 <- entropy(g.one,0.01)
g2 <- entropy(g.one,0.5)
g3 <-entropy(g.one,1.0)

# Simulate a random graph network with n = n of igraph network
n <- vcount(g.one)
attachment <- sum(degree(g.one))/vcount(g.one)
g1.rand <- erdos.renyi.game(n, attachment/(n-1))

g1.rand.o <- entropy(g1.rand,0.01)
g1.rand.t <- entropy(g1.rand,0.5)
g1.rand.th <- entropy(g1.rand,1.0)

# Compare entropy of these random graphs to their original counterparts
par(mfrow=c(1,1))
plot(x=c(0.01,0.5,1.0), y=c(g1,g2,g3), col="red", type="o", pch="o",
     xlab="Beta Values", ylab="Entropy",
     main="Entropies for the igraph Karate network vs
     a related Erdos-Renyi random network",
     ylim=c(min(g1,g2,g3,g1.rand.o,g1.rand.t,g1.rand.th),
          max(g1,g2,g3,g1.rand.o,g1.rand.t,g1.rand.th)), lty=1)
points(x=c(0.01,0.5,1.0), y=c(g1.rand.o,g1.rand.t,g1.rand.th),pch="*",
       col="black")
lines(x=c(0.01,0.5,1.0), y=c(g1.rand.o,g1.rand.t,g1.rand.th), lty=2,
     col="black")
legend("bottomleft",
       legend=c("Karate", "Erdos-Renyi"), col=c("red","black"), pch=c("o","*"),
       lty=c(1,3),ncol=1)