template_funciones_2.py

import numpy as np
import matplotlib.colors as mcolors
import matplotlib.pyplot as plt
import pandas as pd # Para leer archivos
import geopandas as gpd # Para hacer cosas geográficas
import seaborn as sns # Para hacer plots lindos
import networkx as nx # Construcción de la red en NetworkX
import scipy
import template_funciones_1 as tf1

#Matriz A de ejemplo:
    
A_ejemplo = np.array([
    [0, 1, 1, 1, 0, 0, 0, 0],
    [1, 0, 1, 1, 0, 0, 0, 0],
    [1, 1, 0, 1, 0, 1, 0, 0],
    [1, 1, 1, 0, 1, 0, 0, 0],
    [0, 0, 0, 1, 0, 1, 1, 1],
    [0, 0, 1, 0, 1, 0, 1, 1],
    [0, 0, 0, 0, 1, 1, 0, 1],
    [0, 0, 0, 0, 1, 1, 1, 0]
])

s_ejemplo = np.array([1, 1, 1, 1,-1,-1,-1,-1])

def calcula_L(A): # La función recibe la matriz de adyacencia A y calcula la matriz laplaciana
    n = A.shape[1]
    K = np.zeros((n,n))
    for i in range(n): # Calcula la matriz K, que tiene en su diagonal la suma por filas de A
        fila = np.sum(A[i, :])
        K[i][i]=fila
    L = K - A    
    return L

def calcula_R(A): # La función recibe la matriz de adyacencia A y calcula la matriz de modularidad
    E = 0 #número total de conexiones en la red
    n = A.shape[1]
    k = []
    P = np.zeros((n,n))
    for i in range(n): #calculamos los elementos de k y el número total de conexiones E
        fila = np.sum(A[i, :])
        k.append(fila)
        E+=fila    
    for i in range(n): #calculamos P
        for j in range(n):
            P[i][j] = (k[i]*k[j])/(E)  
    R = A - P         
    return R

def calcula_s(v):# Recibe v y retorna s
    n = v.shape[0]
    s = np.zeros(n)
    for i in range(n): #para que cada elemento de s sea el signo de v, dividimos cada elemento de v por su módulo
        s[i]=v[i]/abs(v[i])
    return s 

def calcula_lambda(L,v):# Recibe L y v y retorna el corte asociado
    s = calcula_s(v) 
    lambdon = (1/4)*(s.T)@L@s   
    return lambdon

def calcula_Q(R,v): # La función recibe R y s y retorna la modularidad (a menos de un factor 2E)
    s = calcula_s(v) 
    Q = (s.T)@R@s    
    return Q

def normalizar(v):
   norma2=0
   n=v.shape[0]
   for i in range(n):
       norma2 += v[i]*v[i]
   norma2 = np.sqrt(norma2)
   return v/norma2

def metpot1(A,tol=1e-8,maxrep=np.inf): # Recibe una matriz A y calcula su autovalor de mayor módulo, con un error relativo menor a tol y-o haciendo como mucho maxrep repeticiones
   n = A.shape[0]
   v = np.random.uniform(-1, 1, n) # Generamos un vector de partida aleatorio, entre -1 y 1
   v = normalizar(v) # Lo normalizamos 
   v1 = A@v # Aplicamos la matriz una vez
   v1 = normalizar(v1) # normalizamos
   l = (v.T)@A@v # Calculamos el autovector estimado
   l1 = (v1.T)@A@v1 # Y el estimado en el siguiente paso
   nrep = 0 # Contador
   while np.abs((l1-l))/np.abs(l) > tol and nrep < maxrep: # Si estamos por debajo de la tolerancia buscada 
      v = v1 # actualizamos v y repetimos
      l = l1
      v1 = A@v # Calculo nuevo v1
      v1 = normalizar(v1) # Normalizo
      l = (v.T)@A@v # Calculamos el autovector estimado
      l1 = (v1.T)@A@v1  #Calculo autovector
      nrep += 1 # Un pasito mas
   if not nrep < maxrep:
      print('MaxRep alcanzado')
   l = ((v1.T)@A@v1)/((v1.T)@v1) # Calculamos el autovalor
   return v1,l,nrep<maxrep


def deflaciona(A,tol=1e-8,maxrep=np.inf): # Recibe la matriz A, una tolerancia para el método de la potencia, y un número máximo de repeticiones
    v1,l1,_ = metpot1(A,tol,maxrep) # Buscamos primer autovector con método de la potencia
    deflA = A - (l1 * np.outer(v1, v1)) # Sugerencia, usar la funcion outer de numpy
    return deflA

def metpotI(A,mu,tol=1e-8,maxrep=np.inf): # Retorna el primer autovalor de la inversa de A + mu * I, junto a su autovector y si el método convergió. 
    n = A.shape[0]
    I =np.eye(n,n)
    A = A + (mu*I)
    Ainv = tf1.inversa(A)
    deflA = deflaciona(Ainv)
    return metpot1(deflA,tol=tol,maxrep=maxrep)

def metpotI2(A,mu,tol=1e-8,maxrep=np.inf):
   # Recibe la matriz A, y un valor mu y retorna el segundo autovalor y autovector de la matriz A, 
   # suponiendo que sus autovalores son positivos excepto por el menor que es igual a 0
   # Retorna el segundo autovector, su autovalor, y si el metodo llegó a converger.
   n = A.shape[0]
   I = np.eye(n,n)
   X = A + (mu*I) # Calculamos la matriz A shifteada en mu
   iX = tf1.inversa(X) # La invertimos
   defliX = deflaciona(iX,tol,maxrep) # La deflacionamos
   v,l,_ = metpotI(defliX,mu,tol,maxrep) # Buscamos su segundo autovector
   l = 1/l # Reobtenemos el autovalor correcto
   l -= mu
   return v,l,_


def metpot2(A,v1,l1,tol=1e-8,maxrep=np.Inf):
   # La funcion aplica el metodo de la potencia para buscar el segundo autovalor de A, suponiendo que sus autovectores son ortogonales
   # v1 y l1 son los primeors autovectores y autovalores de A}
   # Have fun!
   return metpot1(deflA,tol,maxrep)

def cortar_matriz(M,v):
    negativos = []
    positivos = []
    s = calcula_s(v)
    n = s.shape[0]
    m=0
    p=0
    nodosp=[]
    nodosm=[]
    for i in range(n): 
         if s[i]>0:
            nodosp.append(i)
            p+=1
            positivos+=[M[i]] #metemos el nodo asociado a s[i] positivo
         else:
            nodosm.append(i)
            m+=1
            negativos+=[M[i]] #metemos el nodo asociado a s[i] negativo
    #convertimos las listas de listas a arrays de numpy:
    positivos=np.array(positivos)
    negativos=np.array(negativos)
    #definimos las dimensiones de Ap y Am:
    Am = np.zeros((m,m))
    Ap = np.zeros((p,p))
    for i in range(p):
        fila=[]
        for n in nodosp:
           fila.append(positivos[i][n]) 
        Ap[i]=fila   
    for i in range(m):
        fila=[]
        for n in nodosm:
           fila.append(negativos[i][n]) 
        Am[i]=fila  
    return Ap,Am 

def laplaciano_iterativo(A,niveles,nombres_s=None):
    # Recibe una matriz A, una cantidad de niveles sobre los que hacer cortes, y los nombres de los nodos
    # Retorna una lista con conjuntos de nodos representando las comunidades.
    # La función debe, recursivamente, ir realizando cortes y reduciendo en 1 el número de niveles hasta llegar a 0 y retornar.
    if nombres_s is None: # Si no se proveyeron nombres, los asignamos poniendo del 0 al N-1
        nombres_s = range(A.shape[0])
    if A.shape[0] == 1 or niveles == 0: # Si llegamos al último paso, retornamos los nombres en una lista
        return([nombres_s])
    else: # Sino:
        L = calcula_L(A) # Recalculamos el L
        v,l,_ = metpotI2(L,1,tol=1e-8,maxrep=np.Inf) # Encontramos el segundo autovector de L
        # Recortamos A en dos partes, la que está asociada a el signo positivo de v y la que está asociada al negativo
        Ap,_ = cortar_matriz(A,v) # Asociado al signo positivo
        _, Am = cortar_matriz(A,v) # Asociado al signo negativo
        return(
                laplaciano_iterativo(Ap,niveles-1,
                                     nombres_s=[ni for ni,vi in zip(nombres_s,v) if vi>0]) +
                laplaciano_iterativo(Am,niveles-1,
                                     nombres_s=[ni for ni,vi in zip(nombres_s,v) if vi<0])
                )        


                         
def modularidad_iterativo(A,R=None,nombres_s=None):
    # Recibe una matriz A, una matriz R de modularidad, y los nombres de los nodos
    # Retorna una lista con conjuntos de nodos representando las comunidades.
    if A is None and R is None:
        print('Dame una matriz')
        return(np.nan)
    if R is None:
        R = calcula_R(A)
    if nombres_s is None:
        nombres_s = range(R.shape[0])
    # Acá empieza lo bueno
    if R.shape[0] == 1: # Si llegamos al último nivel
        return ([nombres_s]) 
    else:
        v,l,_ = metpot1(R,tol=1e-8,maxrep=np.inf) # Primer autovector y autovalor de R
        # Modularidad Actual:
        Q0 = np.sum(R[v>0,:][:,v>0]) + np.sum(R[v<0,:][:,v<0])
        if Q0<=0 or all(v>0) or all(v<0): # Si la modularidad actual es menor a cero, o no se propone una partición, terminamos
            return ([nombres_s])
        else:
            ## Hacemos como con L, pero usando directamente R para poder mantener siempre la misma matriz de modularidad
            Rp,_ = cortar_matriz(R,v) # Parte de R asociada a los valores positivos de v
            _,Rm = cortar_matriz(R,v) # Parte asociada a los valores negativos de v
            vp,lp,_ = metpot1(Rp,tol=1e-8,maxrep=np.inf)  # autovector principal de Rp
            vm,lm,_ = metpot1(Rm,tol=1e-8,maxrep=np.inf) # autovector principal de Rm
        
            # Calculamos el cambio en Q que se produciría al hacer esta partición
            Q1 = 0
            if not all(vp>0) or all(vp<0):
               Q1 = np.sum(Rp[vp>0,:][:,vp>0]) + np.sum(Rp[vp<0,:][:,vp<0])
            if not all(vm>0) or all(vm<0):
                Q1 += np.sum(Rm[vm>0,:][:,vm>0]) + np.sum(Rm[vm<0,:][:,vm<0])
            if Q0 >= Q1: # Si al partir obtuvimos un Q menor, devolvemos la última partición que hicimos
                return([[ni for ni,vi in zip(nombres_s,v) if vi>0],[ni for ni,vi in zip(nombres_s,v) if vi<0]])
            else:
                # Sino, repetimos para los subniveles
                return(modularidad_iterativo(Rp,nombres_s=[ni for ni, vi in zip(nombres_s, v) if vi > 0])+
                       modularidad_iterativo(Rm,nombres_s=[ni for ni, vi in zip(nombres_s, v) if vi < 0]))



def graficar_red2(barrios, museos, m, D):
  A = tf1.construye_adyacencia(D,m) 
  A = np.ceil((1/2)*(A + A.T)) #simetrizamos A
  G = nx.from_numpy_array(A) #construimos la red a partir de la matriz de adyacencia
  modularidad = modularidad_iterativo(A,None)
  n = len(modularidad)
  n=np.log2(n)
  if (n%10>=5):
     n=np.ceil(n)
  else:
     n=np.floor(n)   
  laplaciano = laplaciano_iterativo(A,n,None)

  #construimos el layout con las coordenadas geográficas:
  G_layout = {i:v for i,v in enumerate(zip(museos.to_crs("EPSG:22184").get_coordinates()['x'],museos.to_crs("EPSG:22184").get_coordinates()['y']))}
  factor_escala = 60 #escalamos bastante los nodos para que sean bien visibles
  fig, ax = plt.subplots(1,2, figsize=(8, 8)) #visualiza las redes en el mapa

  comunidades=[laplaciano] + [modularidad] 
  for i,c in enumerate(comunidades): 
    if i==0:
       nombre="laplaciano"
    else:
       nombre="modularidad"
    cant_comunidades = len(c)      
    barrios.to_crs("EPSG:22184").boundary.plot(color='gray',ax=ax[i])
    labels = {n: str(n) if i in c else "" for i, n in enumerate(G.nodes)} #nombres para los nodos
    colores = plt.cm.tab20.colors
    for j,nodo in enumerate(c): #le asignamos un color a cada grupo de nodos
        nx.draw_networkx(G,G_layout,node_size = 2*factor_escala, nodelist = nodo, node_color = [colores[(j%len(colores))]] , ax=ax[i], with_labels=False) # Graficamos red
    nx.draw_networkx_labels(G, G_layout, labels=labels, font_size=10, font_color="k") #agregamos los nombres
    ax[i].set_title(f'{nombre}, m = {m}, {cant_comunidades} comunidades')
  plt.show()