---
title: "Graficas Aleatorias, Redes Sociales y el Internet: Simulaciones"
subtitle: "Mathematics Sin Fronteras, 2021"
author: "Jose A. Sanchez & Mariana Olvera Cravioto"
date: "09/23/2021"
output:
  html_document: default
  pdf_document: default
---

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

## Paquete para visualizar graficas estaticas.
library(igraph)

## Paquete para visualizar graficas interactivamente.
library(networkD3)

## Paquete para visualizar matrices de adyacencia.
library(plot.matrix)

## Otros paquetes y objectos utiles.
if(!require(EnvStats)){
  install.packages("EnvStats")
  library(EnvStats)
}

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



  Generamos distintas visualizaciones de graficas aleatorias a traves de codigo en el lenguaje de programacion R. En particular, generamos visualizaciones para el modelo de Erdos-Renyi, modelo no homogeneo de Erdos-Renyi, modelo en bloques estocastico, grafica de interseccion y conexion preferencial. Para cada grafica, generamos 2 visualizaciones: estatica e interactiva. 
  
  Clarificacion: por cuestiones tecnicas, el texto no tiene tildes. Una disculpa.






## Modelo Erdos-Renyi Homogeneo.

El modelo original de Erdos-Renyi depende de dos parametros: el numero de nodos $n$ y la probabilidad de conexion $p$. Basicamente, incluimos la arista $(i,j)$ en la grafica con probabilidad $p$, y la presencia de las aristas es completamente independiente. Es decir, para $1\leq i<j \leq n$,

$$ \mathbb{P}[(i,j)\in E] = p.$$
Veamos como simular una grafica de Erdos-Renyi en el siguiente bloque de codigo.

```{r, echo = TRUE, fig.height= 6, fig.width= 8, fig.align='center', warning=FALSE}

#########################################################
## Paso 1: generar matriz de adyacencia.
#########################################################

n = 100   # n = total de vertices.
c = 3     # c = grado promedio.
p = c/n   # p = probabilidad de conexion.

adj = matrix(0, n, n)     # Creamos matrix de adyacencia n x n
                          # Por ahora, cero conexiones.

for(.i in 2:n){
  for(.j in 1:(.i-1)){    # Para cada pareja de vertices (.i, .j):
    adj[.i, .j] = rbinom(1,1, p) # Simulamos arista con probabilidad p. 
    
  }
}
adjsymm = adj + t(adj)

#########################################################
## Paso 2: visualizar la grafica con igraph.
#########################################################

ErdosRenyi = graph_from_adjacency_matrix(adjsymm) # Este comando crea un objeto grafica en R.


V(ErdosRenyi)$size <- 5                       # Especificamos algunos parametros visuales 
V(ErdosRenyi)$color <- rep("orange", n)       # para nuestra grafica:
V(ErdosRenyi)$label <- ""                     
V(ErdosRenyi)$frame.color <- "white"          
E(ErdosRenyi)$arrow.mode <- 0
lw <- layout_with_fr(ErdosRenyi)

title = paste0("Erdos-Renyi Random Graph ",   # Titulo de la visualizacion.
               "(n = ", n, 
               ", p = ", p, ")")

plot(ErdosRenyi, layout = lw, main = title)


#########################################################
## Paso 3: Visualizacion interactiva.
#########################################################

edges = which(adj == 1, arr.ind = TRUE)
edges = as.data.frame(edges)

p <- simpleNetwork(edges, height="100px", width="100px",        
        Source = 1,                 
        Target = 2,                 
        linkDistance = 10,          
        charge = -100,                
        fontSize = 1,               
        fontFamily = "serif",       
        linkColour = "#666",        
        nodeColour = "#69b3a2",     
        opacity = 0.9,              
        zoom = T                    
        )
p

#########################################################
## Paso 4: Juega con los parametros
#########################################################

###### n = 500, p = 2/n
## Esperamos una componente conexa gigante.

###### n = 500, p = 1/n
## Esperamos una componente de tama�o 500^(2/3) = 62.

###### n = 500, p = 0.5/n
## Esperamos una grafica casi completamente disconexa.

```





\newpage



## Modelo de Erdos-Renyi No-Homogeneo

En el caso de un modelo no homogeneo, cada vertice $i$ tiene un peso asociado $w_i$. Entre mas grande sea el peso de un vertice, mas popular el vertice sera. Por lo tanto, tendra mas conexiones. Matematicamente, si $L= \sum_{i=1}^n w_i$, entonces,
$$ \mathbb{P}[(i,j)\in E] = \min\left\{\frac{w_iw_j}{L}, 1\right\}$$
En este caso, el grado promedio del vertice $i$ es aproximadamente $w_i$.

```{r, fig.width=8, fig.height=6, fig.align='center', warning=FALSE}

#########################################################
## Paso 1: generar matriz de adyacencia.
#########################################################

n = 100                     # Numero de vertices. 


b = 1.5                       # Simulacion de la pareto.
alpha = 2
U = runif(n = n, min = 0, max =  1)
w = b/( (1 - U)^(1/alpha) ) # Pesos!
L = sum(w)                  # Suma de todos los pesos.


adj = matrix(0, n, n)       # Creamos matrix de adyacencia n x n
                            # Por ahora, cero conexiones.

for(.i in 2:n){                         
  for(.j in 1:(.i-1)){          # Para cada pareja de vertices (.i, .j):
    .p = min(w[.i]*w[.j]/L, 1)    # a) Encontramos probabilidad de conexion.
    adj[.i, .j] = rbinom(1,1, .p) # b) Simulamos arista con probabilidad p. 
  }
}
adjsymm = adj + t(adj)

#########################################################
## Paso 2: visualizar la grafica con igraph.
#########################################################


NH.ErdosRenyi = graph_from_adjacency_matrix(adjsymm)  # Creamos grafica con igraph.
V(NH.ErdosRenyi)$size <- 6
V(NH.ErdosRenyi)$label <- ""                       
V(NH.ErdosRenyi)$frame.color <- "white"
E(NH.ErdosRenyi)$arrow.mode <- 0


## Quiero colorear vertices de acuerdo a sus grados.
deg = apply(adjsymm, MARGIN = 1, sum)         # a) Calculamos los grados.
dmax = max(deg) + 1                       # b) Encontramos el grado maximo + 1

colfunc <- colorRampPalette(c("orange",   # c) Generamos un gradiente de color
                              "red",      #     del naranja a negro.
                              "black"))
colors = colfunc(dmax)                    

V(NH.ErdosRenyi)$color <- colors[deg + 1]     # Coloreamos los vertices de acuerdo
                                              # a su grado.

lw <- layout_with_fr(NH.ErdosRenyi)
plot(NH.ErdosRenyi, layout = lw, main = "Inhomogeneous Erdos-Renyi: Pareto weights")



#########################################################
## Paso 3: Visualizacion dinamica
#########################################################

edges = which(adj == 1, arr.ind = TRUE)
edges = as.data.frame(edges)

p <- simpleNetwork(edges, height="100px", width="100px",        
        Source = 1,                 
        Target = 2,                 
        linkDistance = 30,          
        charge = -100,                
        fontSize = 1,               
        fontFamily = "serif",       
        linkColour = "#666",        
        nodeColour = "#69b3a2",     
        opacity = 0.9,              
        zoom = T                    
        )
p




#########################################################
## Paso 3: Juega con los parametros
#########################################################


```





\newpage



## Modelo en Bloques Estocastico

Para el modelo en bloques estocastico, los vertices de la grafica pertenecen a una de $K$ comunidades diferentes.  La probabilidad de que un vertice $i$ en la comunidad $A$ y otro vertice $j$ en la comunidad $B$ se conected esta dada por  $p_{A,B}= \frac{\kappa(A,B)}{n}$, donde $\kappa:\{1,\ldots, K\}^2\to \mathbb{R}_{+}$ es un kernel de conexion.


```{r, fig.width=9, fig.height=6, fig.align='center', warning=FALSE}

#########################################################
## Paso 1: generar matriz de adyacencia.
#########################################################

n1 = 100                 # Creamos 3 comunidades del mismo tama�o.
n2 = 100
n3 = 100
n = n1 + n2 + n3
class = rep(c(1,2,3),    
            c(n1,n2,n3))

cmat = matrix(c(15, 0.2, 0.2,  # Generamos el kernel.
                0.2, 15, 0.2,
                0.2, 0.2, 15), 
              ncol= 3, 
              byrow = TRUE)
pmat = cmat / n        # Matriz de probabilidad de conexiones entre comunidades.



adj = matrix(0, n, n)       # Creamos matrix de adyacencia n x n
                            # Por ahora, cero conexiones.
for(.i in 2:n){                         
  for(.j in 1:(.i-1)){            # Para cada pareja de vertices (.i, .j):
    .p  = pmat[class[.i], class[.j]]  # a) Encontramos la probabilidad de conexion entre clases.
    adj[.i, .j] = rbinom(1,1, .p)     # b) Simulamos arista con probabilidad p. 
  }
}
adjsymm = adj + t(adj)


#########################################################
## Paso 2: visualizar la grafica con igraph.
#########################################################

SBM = graph_from_adjacency_matrix(adjsymm)    # Crear grafica con i-graph.
V(SBM)$size <- 6
V(SBM)$label <- "" 
V(SBM)$frame.color <- "white"
E(SBM)$arrow.mode <- 0

## Para colorear cada vertice de acuerdo a su pertenencia.
V(SBM)$color <- rep(cbPalette[c(2,3,4)], c(n1,n2,n3))


lw <- layout_with_fr(SBM)

plot(SBM, layout = lw, main = "Stochastic Block Model: 3 communities")



#########################################################
## Paso 3: Grafica interactiva.
#########################################################

edges = which(adj == 1, arr.ind = TRUE)
edges = as.data.frame(edges)

p <- simpleNetwork(edges, height="100px", width="100px",        
        Source = 1,                 
        Target = 2,                 
        linkDistance = 10,          
        charge = -100,                
        fontSize = 1,               
        fontFamily = "serif",       
        linkColour = "#666",        
        nodeColour = "#69b3a2",     
        opacity = 0.9,              
        zoom = T                    
        )
p





#########################################################
## Paso 3: Juega con los parametros
#########################################################

####### Cosas a modificar:

## La matriz de conexiones por comunidades pmat.
## 1) �Que pasa si incrementamos la probabilidad entre comunidades?
## 2) �Que pasa si damos diferentes probabilidades de conexion para las comunidades.

## Las comunidades.
## 3) �Que pasa si las comunidades se vuelven mas grandes?
## 4) �Que pasa si las comunidades no son del mismo tama�o?


```



\newpage





### Grafica con Clusters.

El resto de las graficas que hemos visualizado hasta ahora tienden a tener "ciclos largos." Con esto, nos referimos al hecho de que, la longitud de un ciclo en promedio es $O(\log(n))$, donde $n$ es el numero de vertices. Esto significa que, conforme $n$ crece, los ciclos en la grafica se vuelven progresivamente mas largos, y la presencia de triangulos se vuelve cada vez menos probable.

Esto puede ser poco realista, puesto que en la vida real, es comun que haya muchos triangulos en redes. Los triangulos corresponden a "amistades en comun" lo cual ocurre constantemente en redes sociales. 

El modelo de interseccion genera graficas aleatorias con una alta presencia de triangulos.

```{r, fig.width= 10, fig.height=8, fig.align='center', warning=FALSE}

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

n = 100                     # Numero de vertices (actores).
beta = 0.5
gamma = 0.5
m = floor(beta*n)          # Numero de vertices opuestos (peliculas)

# old/good.
#b = 3                       # Simulacion de la pareto.
#alpha = 2.5

b = 5                       # Simulacion de la pareto.
alpha = 3.5
U = runif(n = n, min = 0, max =  1)
w = b/( (1 - U)^(1/alpha) ) # Pesos!
L = sum(w)                  # Suma de todos los pesos.
  

bip = matrix(0, ncol = n2, nrow = n)
for(.i in 1:n){             
  for(.j in 1:m){                  # Para cada par (actor, pelicula):
    .pi = min(gamma*w[.i]/n, 1)       # a) Encontramos probabilidad de conexion.
    bip[.i, .j] = rbinom(n = 1,       # b) Generamos arista con probabilidad pi
                         size = 1, 
                         prob = .pi)
  }
}

#########################################################
## Paso 2: generar la grafica de peliculas en comun.
#########################################################

#edges = matrix(0, nrow = 1, ncol = 2)
adj = matrix(0, n, n)
count = 0
for(.i in 2:n){
  for(.j in 1:(.i-1)){          # Para cada par de actores:
    
    .iftrian = sum(bip[.i, ] * bip[.j, ]) # iftrian = "�Cuantas cuantas peliculas
                                          #             han colaborado i y j?"
    
    if(.iftrian > 0){   # Si han colaborado, conectar, asociar en la grafica.
      adj[.i, .j] = 1
    } 
  } 
}
adjsymm = adj + t(adj)

#########################################################
## Paso 3: generar la grafica de peliculas en comun.
#########################################################

IntGraph = graph_from_adjacency_matrix(adjsymm)
V(IntGraph)$size <- 4
V(IntGraph)$label <- "" 
V(IntGraph)$frame.color <- "white"
E(IntGraph)$arrow.mode <- 0

## Quiero colorear vertices de acuerdo a sus grados.
{
  deg = apply(adjsymm, MARGIN = 1, sum)         # a) Calculamos los grados.
  dmax = max(deg) + 1                       # b) Encontramos el grado maximo + 1
  
  colfunc <- colorRampPalette(c("orange",   # c) Generamos un gradiente de color
                                "red",      #     del naranja a negro.
                                "black"))
  colors = colfunc(dmax)                    
  
  V(IntGraph)$color <- colors[deg + 1]     # Coloreamos los vertices de acuerdo
                                                # a su grado.
}

lw <- layout_with_fr(IntGraph)
coeff = transitivity(IntGraph, type = "average")
plot(IntGraph, layout = lw, 
     main = paste("Intersection model: clustering = ", coeff) )



#########################################################
## Paso 4: Visualizacion dinamica
#########################################################

edges = which(adj == 1, arr.ind = TRUE)
edges = as.data.frame(edges)

p <- simpleNetwork(edges, height="100px", width="100px",        
        Source = 1,                 
        Target = 2,                 
        linkDistance = 30,          
        charge = -100,                
        fontSize = 1,               
        fontFamily = "serif",       
        linkColour = "#666",        
        nodeColour = "#69b3a2",     
        opacity = 0.9,              
        zoom = T                    
        )
p



```



\newpage





## Conexion preferencial.

En el modelo de conexion preferencial, los vertices llegan a la grafica en secuencia, y se conectan a los vertices con mas conexiones con mayor probabilidad. Por esto, los vertices mas viejos (que llegan primero) tienden a tener muchas conexiones, comparados con los ultimos.

```{r, fig.width= 10, fig.height=8, fig.align='center', warning=FALSE}

#########################################################
## Paso 1: generar la matriz de adyacencia.
#########################################################

n = 100
deg = rep(c(2,1), c(1,n-1))

adj = diag(n)
adjsymm = diag(n)

for(.k in 2:n){
  .ps = deg[1:.k] / sum(deg[1:.k])
  adj[.k, .k] = 0
  adjsymm[.k, .k] = 0
  
  .j = sample(x = 1:.k, size = 1, prob = .ps) 
  
  adj[.k, .j] = 1
  adjsymm[.k, .j] = 1
  adjsymm[.j, .k] = 1
    
  deg = apply(adjsymm, MARGIN = 1, sum)
}


#########################################################
## Paso 3: generar la grafica de peliculas en comun.
#########################################################

AlbBar = graph_from_adjacency_matrix(adjsymm)
V(AlbBar)$size <- 4 
V(AlbBar)$frame.color <- "white"
E(AlbBar)$arrow.mode <- 0

lw <- layout_with_fr(AlbBar)
coeff = transitivity(AlbBar, type = "average")
plot(AlbBar, layout = lw, 
     main = paste("Preferential Attachment"))



#########################################################
## Paso 4: Visualizacion interactiva.
#########################################################

edges = which(adj == 1, arr.ind = TRUE)
edges = as.data.frame(edges)

p <- simpleNetwork(edges, height="100px", width="100px",        
        Source = 1,                 
        Target = 2,                 
        linkDistance = 40,          
        charge = -50,                
        fontSize = 15,               
        fontFamily = "serif",       
        linkColour = "#666",        
        nodeColour = "#69b3a2",     
        opacity = 0.9,              
        zoom = T                    
        )
p


```









