####SOME EXPLORATIVE NETWORK FUNCTIONS IN STATNET AND IGRAPH
############################################################################ 

library(statnet)
library(igraph)

load("netgot_data.RData")

## Global Network indexes (statnet)

netgot
network.size(netgot)
isolates(netgot)


netgot_friend
network.size(netgot_friend)
isolates(netgot_friend)
length(isolates(netgot_friend))

netgot_friend <- network::delete.vertices(netgot_friend, isolates(netgot_friend))
network.size(netgot_friend)

netgot_enemy
netgot_enemy<- network::delete.vertices(netgot_enemy, isolates(netgot_enemy))
network.size(netgot_enemy)

gden(netgot, mode="graph", ignore.eval==FALSE) # Density
grecip(netgot, measure="dyadic.nonnull") # Dyadic reciprocity

gtrans(netgot, mode="graph") # Transitivity

gtrans(netgot[,], mode="graph") # Transitivity

gden(netgot_friend, mode="graph") # Density
gden(netgot_enemy, mode="graph") # Density

gtrans(netgot_friend[,], mode="graph") # Transitivity
gtrans(netgot_enemy[,], mode="graph") # Transitivity


############################################################################ 
## centrality: degree, closeness, betweenness
               
# degree centrality
               
deg <- sna::degree(netgot, gmode="graph") # Default: total degree
ideg <- sna::degree(netgot, cmode="indegree") 
odeg <- sna::degree(netgot, cmode="outdegree")

?degree
# Once centrality scores are computed, we can handle them using standard R methods:

plot(ideg, odeg, type="n", xlab="inDegree (popularity)", ylab="outDegree (activity)") # Plot ideg by odeg
abline(0, 1, lty=3)
text(jitter(ideg), jitter(odeg), network.vertex.names(netgot), cex=0.75, col=2)

max(deg)

# Let see what are the degree in the friend and enemy networks 

degF <- sna::degree(netgot_friend) # total friend degrees
 
degE <- sna::degree(netgot_enemy) # total enemy degrees


#Plot simple histograms of the degree distribution:
par(mfrow=c(1,2)) # Set up a 2x2 display
hist(degF, xlab="degree", main="Friends degree distribution", prob=TRUE)
hist(degE, xlab="degree", main="Enemy degree distribution", prob=TRUE)
par(mfrow=c(1,1)) # Restore display

dev.off()
# 
plot(degF, degE, type="n", xlab="Degree (friends)", ylab="Degree (enemy)") 
abline(0, 1, lty=3)
text(jitter(degF), jitter(degE), network.vertex.names(netgot), cex=0.75, col=2)

               
# Centrality scores can also be used with other sna routines, e.g., \texttt{gplot()}:
pdf("degreeGoT.pdf", width = 35, height = 35)
gplot(netgot, vertex.cex=deg^0.5/2, vertex.sides=1000,label.cex=.5, gmode="graph", displaylabels=TRUE, displayisolates = F)
dev.off()

gplot(netgot_friend, vertex.cex=degF^0.5/2, vertex.sides=1000, label.cex=.5, gmode="graph", displaylabels=TRUE, displayisolates = F)

gplot(netgot)

#Betweenness & closeness:

?betweenness
bet <- sna::betweenness(netgot, gmode="graph") # betweenness
bet

betdf <- data.frame(netgot %v% "name", bet) 
betdf[order(betdf$bet, decreasing = TRUE), c(1,2)]

gplot(netgot, vertex.cex=sqrt(bet)/10, gmode="graph", displaylabels = T, label.cex = 0.8) # Use with gplot

clo <- sna::closeness(netgot) # closeness

hist(bet)

?evcent

ev <- sna::evcent(netgot) # eigenvector centrality
ev 

evdf <- data.frame(netgot %v% "name",ev) 

evdf[order(evdf$ev, decreasing = TRUE), c(1,2)]

?centralization
centralization(netgot, degree) # Concentration of edges
           
centralization(netgot, mode="graph", sna::closeness)

centralization(netgot, sna::betweenness)

    
############################################################################ 
##  Subnetwork censuses

sna::dyad.census(netgot) # M,A,N counts
sna::dyad.census(netgot_friend)
sna::dyad.census(netgot_enemy)

## triad census

sna::triad.census(netgot)
triad.census(netgot_friend)
triad.census(netgot_enemy)

triad.census(netgot, mode="graph") # Undirected triad census

tc<-sna::triad.census(netgot) # directed triad census
barplot(tc, names.arg=colnames(tc), las=2)

barplot(tc[-1]/sum(triad.census(netgot)[-1]), names.arg=colnames(tc)[-1], las=2)

#Component information can be obtained in various ways
               
comp <- sna::components(netgot) # Strong component count
comp <- components(netgot_friend, connected="weak") # Weak component count

gplot(netgot_friend)

cd <- component.dist(netgot_friend, connected="weak") # Get weak components
max(cd$csize)/network.size(netgot_friend)

# Who's in the largest component?
               
cl <- component.largest(netgot_friend, connected="weak")
cl
               
# Plot the largest component

gplot(netgot_friend[cl,cl], boxed.lab=FALSE, label.cex=0.5,label.col=4,label=network.vertex.names(netgot)[cl])

#cliques census

cc<-clique.census(netgot_friend,clique.comembership="bysize")

############################################################################ 
## A Smaller net example with igraph

install.packages("sand")
library(sand)

data("karate")
karate

vcount(karate)
ecount(karate)

edge_density(karate)


plot(karate)

hist(degree(karate), col="lightblue", xlim=c(0, 50), xlab="Vertex Degree", ylab="Frequency", main="") #degree

hist(graph.strength(karate), col="pink", xlab="Vertex Strength", ylab="Frequency") #weigted deg

d.karate <- degree(karate)

dd.karate <- degree.distribution(karate)
?degree.distribution

############################
#### Edge scores (betweenness)

eb <- edge.betweenness(karate)
E(karate)[order(eb, decreasing=T)[1:5]]

############################
#### Clique census

cliques.karate <- cliques(karate)

which(sapply(cliques(karate), length)==4)

induced.subgraph(karate, cliques.karate[[159]])

plot(induced.subgraph(karate, cliques.k[[103]]))

table(sapply(cliques(karate), length))

cliques(karate)[sapply(cliques(karate), length) == 5] ### note the redundancy

# Cliques of larger sizes necessarily include cliques of smaller sizes. 
# A maximal clique is a clique that is not a subset of a larger clique. 
# In the karate network, only two of the 11 cliques of size four are maximal

table(sapply(maximal.cliques(karate), length))

############################
#### Ego-density (or local density)


ego.instr <- induced.subgraph(karate, neighborhood(karate, order= 1, nodes= "Mr Hi")[[1]])
plot(ego.instr)

ego.admin <- induced.subgraph(karate, neighborhood(karate, order= 1, nodes= "John A")[[1]])


edge_density(karate)

edge_density(ego.instr)

edge_density(ego.admin)

############################
### Other centrality

closeness(karate)
betweenness(karate)
betweenness(karate, normalized=T)
evcent(karate)$vector

?betweenness

############################
### Transitivity

transitivity(karate)

transitivity(karate, "local", vids=c(1,34))
transitivity(karate, "local", vids=c("Mr Hi", "John A"))


###########################
#### path length

average.path.length(karate)
diameter(karate)

###########################
#### components

comps <- decompose.graph(karate)
table(sapply(comps, vcount))

karate.gc <- decompose.graph(karate)[[1]]
       plot(karate.gc)
       
data(aidsblog)
?aidsblog
plot(aidsblog, vertex.label="", vertex.size=5)

data(yeast)
plot(yeast, vertex.label="", vertex.size=5)
?yeast

vcount(yeast)
ecount(yeast)

edge_density(yeast)

hist(degree(yeast))

comps <- decompose.graph(yeast)
table(sapply(comps, vcount))

yeast.gc <- decompose.graph(yeast)[[1]]
plot(yeast.gc, vertex.label="", vertex.size=5)