---
title: "Random Graphs, Social Networks and the Internet: Centrality"
subtitle: "Mathematics Sin Fronteras, 2021"
author: "Jose A. Sanchez & Mariana Olvera Cravioto"
date: "10/14/2021"
output:
  html_document: default
  pdf_document: default
---

In the following, we generate visualizations of different random graph models. Furthermore, on each graph, we compare different measurements of centrality for vertices: degree centrality, closeness centrality, betweenness centrality and page rank. 


```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

## Package used to save and visualize networks.
library(igraph)

## Package used to visualize networks dynamically.
library(networkD3)

## Package for visualizing adjacency matrices.
library(plot.matrix)

## Not 100% sure of what this is for.
if(!require(EnvStats)){
  install.packages("EnvStats")
  library(EnvStats)
}

cbPalette <- c("#999999", "#E69F00", "#56B4E9", 
               "#009E73", "#F0E442", "#0072B2", 
               "#D55E00", "#CC79A7")


#########################################################
## Paso 2: Measuring centrality:
#########################################################

CentralityNames = c("Degree Centrality", "Closeness Centrality", "Betweenness Centrality", "Page-Rank")

###########################
Laplacian = function(Theta){
  
  .diag = apply(Theta, MARGIN = 1, sum)
  
  .lap = diag(.diag) - Theta
}

###########################
MakeConnected = function(adjsymm, n){
  .Lap = Laplacian(Theta = adjsymm)
  
  .eig = eigen(.Lap)
  .Ncomm = sum(abs(.eig$values) < 0.001)
  
  if(.Ncomm == 1){
    return(adjsymm)
  }
  .adjsymm = adjsymm
  
  .comm = kmeans(x = .eig$vectors[, -(1:(n-.Ncomm))], centers = .Ncomm)
  .comm = .comm$cluster
  
  .main = order(table(.comm), decreasing = TRUE)[1]
  
  for(.clust1 in (1:.Ncomm)[-.main]){
    .clust2 = .main
    
    .members1 = which(.comm == .clust1)
    .members2 = which(.comm == .clust2)
    .nmemb1 = length(.members1)
    .nmemb2 = length(.members2)
    
    
    .rind1 = sample(x = 1:.nmemb1, size = 1)
    .rind2 = sample(x = 1:.nmemb2, size = 1) 
    
    .node1 = .members1[.rind1]
    .node2 = .members2[.rind2]
    
    .adjsymm[.node1, .node2] = 1
    .adjsymm[.node2, .node1] = 1
  }
  return(.adjsymm)
}

###########################
CentralityMeasures = function(graph, n, which = 4){
  .cent = NULL
  if( which == 1 ){
    .cent = degree(graph = graph)
  } else if( which == 2 ){
    .cent = closeness(graph = graph)
  } else if( which == 3 ){
    .cent = betweenness(graph = graph)
  } else if( which == 4 ){
    .cent = page.rank(graph = graph)$vector
  }
  
  return(.cent)
}

```


## Erdos-Renyi Model




```{r pressure, echo=FALSE, fig.width=10, fig.height = 8}

set.seed(2)

#########################################################
## Step 1: create the adjacency matrix.
#########################################################

n = 100   # n = total node count.
c = 3     # c = average degree.
p = c/n   # p = connection probability.

adj = matrix(0, n, n)     # We create an adjacency matrix.
                          # We start with no connections.

for(.i in 2:n){
  for(.j in 1:(.i-1)){    # For each pair of vertices (.i, .j):
    adj[.i, .j] = rbinom(1,1, p) # we add the edge with probability p´.
    
  }
}
adjsymm = adj + t(adj)
adjsymm = MakeConnected(adjsymm = adjsymm, n = n)

graph = graph_from_adjacency_matrix(adjsymm) 


graph = as.undirected(graph, mode = "collapse") # Specify visual parameters.
                  
V(graph)$label <- ""                     
V(graph)$frame.color <- "white"          
E(graph)$arrow.mode <- 0
lw <- layout_with_fr(graph)


#########################################################
## Paso 2: Generate plots of centrality measurements.
#########################################################


for(.which in 1:4){
  ## Find measure of centrality:
  .cent = CentralityMeasures(graph = graph, n = n, which = .which)
  .maxcent = max(.cent)
  .centcol = 100 * .cent/.maxcent
  
  ## Color the vertices of the graph.
  colfunc <- colorRampPalette(c("orange", "red", "black"))
  colors = colfunc(101)

  V(graph)$color <- colors[.centcol + 1]
  V(graph)$size <- sqrt(.centcol) + 1
  
  plot(graph, 
       layout = lw, 
       vertex.label.dist=1,
       main = CentralityNames[.which])

}



```


## Chung-Lu Model


```{r echo=FALSE, fig.width=10, fig.height = 8}

set.seed(4)

#########################################################
## Paso 1: generate graph.
#########################################################

n = 100                     # n = total vertex count.


b = 1.5                       # Generating pareto weights.
alpha = 1.5
U = runif(n = n, min = 0, max =  1)
w = b/( (1 - U)^(1/alpha) ) # Weights!
L = sum(w)                  # Sum of all weights.


adj = matrix(0, n, n)       # We create an n x n adjacency matrix.
                            # We start with no connections.

for(.i in 2:n){                         
  for(.j in 1:(.i-1)){          # For each pair of vertices (.i, .j):
    .p = min(w[.i]*w[.j]/L, 1)    # a) We find the connection probability.
    adj[.i, .j] = rbinom(1,1, .p) # b) We include the edge in the network with probability p. 
  }
}
adjsymm = adj + t(adj)
adjsymm = MakeConnected(adjsymm = adjsymm, n = n)

graph = graph_from_adjacency_matrix(adjsymm) # This command creates the graph object based on the adjacency matrix.
graph = as.undirected(graph, mode = "collapse") # Specify visual parameters.
V(graph)$label <- ""                     
V(graph)$frame.color <- "white"          
E(graph)$arrow.mode <- 0
lw <- layout_with_fr(graph)

#########################################################
## Paso 2: Generate plots of centrality measurements.
#########################################################


for(.which in 1:4){
  
  ## Find measure of centrality:
  .cent = CentralityMeasures(graph = graph, n = n, which = .which)
  .maxcent = max(.cent)
  .centcol = 100 * .cent/.maxcent
  
  ## Color the vertices of the graph.
  colfunc <- colorRampPalette(c("orange", "red", "black"))
  colors = colfunc(101)
  
  V(graph)$color <- colors[.centcol + 1]
  V(graph)$size <- sqrt(.centcol) + 1    
  
  plot(graph, 
       layout = lw, 
       vertex.label.dist=1,
       main = CentralityNames[.which])

}

```




## Intersection Model

```{r echo=FALSE, fig.width=10, fig.height = 8}

set.seed(6)

#########################################################
## Paso 1: generar la grafica bipartita.
#########################################################

n = 100                     # Total vertex count (actores).
beta = 0.5
gamma = 0.5
n2 = floor(beta*n)          # Number of vertices on the bipartite opposite side (movies)

b = 5                       # We simulate Pareto weights.
alpha = 3.5
U = runif(n = n, min = 0, max =  1)
w = b/( (1 - U)^(1/alpha) ) # Weights!!
L = sum(w)                  # Sum of weights.
  

bip = matrix(0, ncol = n2, nrow = n)
for(.i in 1:n){             
  for(.j in 1:n2){                  # for each pair of (actor, movie):
    .pi = min(gamma*w[.i]/n, 1)       # a) we find the connection probability.
    bip[.i, .j] = rbinom(n = 1,       # b) We add a connection between them with probability p.
                         size = 1, 
                         prob = .pi)
  }
}

#########################################################
## Paso 2: Generate the adjacency matrix of actors.
#########################################################


adj = matrix(0, n, n)
count = 0
for(.i in 2:n){
  for(.j in 1:(.i-1)){          # For each pair of actors.
    
    .iftrian = sum(bip[.i, ] * bip[.j, ]) # iftrian = "How many movies did i and j
                                          #             collaborate on?"
    
    if(.iftrian > 0){   #  If i and j have collaborated, connect them.
      adj[.i, .j] = 1
    } 
  } 
}
adjsymm = adj + t(adj)
adjsymm = MakeConnected(adjsymm = adjsymm, n = n)

graph = graph_from_adjacency_matrix(adjsymm)  # This command creates the graph object 
                                              # based on the adjacency matrix.
V(graph)$label <- ""                     
V(graph)$frame.color <- "white"          
E(graph)$arrow.mode <- 0
lw <- layout_with_fr(graph)

#########################################################
## Paso 2: Generate plots of centrality measurements.
#########################################################


for(.which in 1:4){
  
  ## Find measure of centrality:
  .cent = CentralityMeasures(graph = graph, n = n, which = .which)
  .maxcent = max(.cent)
  
  .centcol = 100 * .cent/.maxcent
  
  #V(graph)$label <- round( 100 * .cent/.maxcent , digits = 2)
  
  ## Color the vertices of the graph.
  colfunc <- colorRampPalette(c("orange", "red", "black"))
  colors = colfunc(101)

  V(graph)$color <- colors[.centcol + 1]
  V(graph)$size <- sqrt(.centcol) + 1   
  
  plot(graph, 
       layout = lw, 
       vertex.label.dist=1,
       main = CentralityNames[.which])

}

```