redes - python ping script

¿Cómo calcular Eb(k) de redes con Python?

Question

En el documento titulado Escalado de las correlaciones de grado y su influencia en la difusión en redes sin escala , los autores definen la cantidad de $ E_b (k) $ para medir el grado de correlación de grado.

Papel

LK Gallos, C. Song y HA Makse, escalado de correlaciones de grado y su influencia en la difusión en redes libres de escala, Phys. Rev. Lett. 100, 248701 (2008).

Puede leer el artículo siguiendo este enlace o leer el libro de Google relacionado.

Pregunta

Mi pregunta es cómo calcular Eb (k) de redes con Python? Mi problema es que no puedo reproducir los resultados de los autores. Lo pruebo utilizando los datos de Condense Matter. El resultado de Eb (k) se muestra en la figura de arriba. ¡Puedes ver que un problema en mi figura es que Eb (k) es mucho más grande que 1! !! También he probado Internet (como datos de nivel) y los datos de WWW, y el problema persiste. Sin duda, hay algo seriamente mal con mi algoritmo o código. Puedes reproducir mis resultados y compararlos con los de los autores. Su solución o sugerencia son muy apreciadas. Presentaré mi algoritmo y el script de python a continuación.

Sigo los siguientes pasos:

1. Para cada borde, para encontrar los bordes cuyos k = k, y k ''> 3k. La probabilidad de estos bordes se denota como P (k, k '')
2. Para nodo, para obtener la proporción de nodos cuyo grado es mayor que b * k, que se denota como p (k ''), también podemos tener k'' * p (k '')
3. Para obtener el numerador P1: p1 = / sum P (k, k '') / k'' * P (k '')
4. Para obtener el denominador p2: P2 = / sum P (k '')
5. Eb (k) = p1 / p2

Script de Python

La secuencia de comandos de Python se da a continuación:

%matplotlib inline
import networkx as nx
import matplotlib.cm as cm
import matplotlib.pyplot as plt
from collections import defaultdict

def ebks(g, b):
    edge_dict = defaultdict(lambda: defaultdict(int))
    degree_dict = defaultdict(int)
    edge_degree = [sorted(g.degree(e).values()) for e in g.edges()]
    for e in edge_degree:
        edge_dict[e[0]][e[-1]] +=1
    for i in g.degree().values():
        degree_dict[i] +=1
    edge_number = g.number_of_edges()
    node_number = g.number_of_nodes()
    ebks, ks = [], []
    for k1 in edge_dict:
        p1, p2 = 0, 0
        for k2 in edge_dict[k1]:
            if k2 >= b*k1:
                pkk = float(edge_dict[k1][k2])/edge_number
                pk2 = float(degree_dict[k2])/node_number
                k2pk2 = k2*pk2
                p1 += pkk/k2pk2
        for k in degree_dict:
            if k>=b*k1:
                pk = float(degree_dict[k])/node_number
                p2 += pk
        if p2 > 0:
            ebks.append(p1/p2)
            ks.append(k1)
    return ebks, ks

Pruebo con los datos de ca-CondMat, puede descargarlos de esta url: http://snap.stanford.edu/data/ca-CondMat.html

# Load the data
# Remember to change the file path to your own
ca = nx.Graph()
with open (''/path-of-your-file/ca-CondMat.txt'') as f:
    for line in f:
        if line[0] != ''#'':
            x, y = line.strip().split(''/t'')
            ca.add_edge(x,y)
nx.info(ca)

#calculate ebk 
ebk, k = ebks(ca, b=3)

plt.plot(k,ebk,''r^'')
plt.xlabel(r''$k$'', fontsize = 16)
plt.ylabel(r''$E_b(k)$'', fontsize = 16)
plt.xscale(''log'')
plt.yscale(''log'')
plt.show()

Actualización : El problema aún no se ha resuelto.

def ebkss(g, b, x):
    edge_dict = defaultdict(lambda: defaultdict(int))
    degree_dict = defaultdict(int)
    edge_degree = [sorted(g.degree(e).values()) for e in g.edges()]
    for e in edge_degree:
        edge_dict[e[0]][e[-1]] +=1
    for i in g.degree().values():
        degree_dict[i] +=1
    edge_number = g.number_of_edges()
    node_number = g.number_of_nodes()
    ebks, ks = [], []
    for k1 in edge_dict:
        p1, p2 = 0, 0
        nk2k = np.sum(edge_dict[k1].values())
        pk1 = float(degree_dict[k1])/node_number
        k1pk1 = k1*pk1
        for k2 in edge_dict[k1]:
            if k2 >= b*k1:
                pk2k = float(edge_dict[k1][k2])/nk2k
                pk2 = float(degree_dict[k2])/node_number
                k2pk2 = k2*pk2
                p1 += (pk2k*k1pk1)/k2pk2
        for k in degree_dict:
            if k>=b*k1:
                pk = float(degree_dict[k])/node_number
                p2 += pk
        if p2 > 0:
            ebks.append(p1/p2**x)
            ks.append(k1)
    return ebks, ks

Answer 1

Parece que en realidad estás calculando la probabilidad condicional usando distribuciones discretas, por lo que estás obteniendo muchos ceros, lo que crea problemas.

En el documento (parte superior de la segunda columna, segunda página), parece que están usando una ley de potencia ajustada a los datos para reemplazar los valores discretos ruidosos con una función suave y agradable. Y eso es también, supongo, por qué escriben E_b en términos de integrales en lugar de sumas.

Si yo fuera usted, les pediría su código a los autores del artículo. Y luego le pediría a la revista que deje de publicar artículos sin el código de respaldo.

Answer 2

Según el documento, el propósito de Eb (k) es obtener el exponente de correlación épsilon: "[Introducimos] una cantidad invariante de escala Ebk para simplificar la estimación de épsilon" (segunda página, parte inferior de la primera columna).

No he encontrado una manera de hacer que Eb (k) <1, pero he encontrado una corrección que calcula correctamente épsilon .

De acuerdo con la ecuación 4, Eb (k) ~ k ^ - (epsilon-gamma) (donde la distribución de grados P (k) ~ k ^ -gamma, una ley de poder). Por lo tanto, si trazamos la pendiente de log (Eb (k)) contra log (k), deberíamos obtener gamma - epsilon. Sabiendo gamma, podemos fácilmente obtener épsilon.

Tenga en cuenta que esta pendiente es invariable si Eb (k) está escalada por una constante. Por lo tanto, el problema con su Eb (k) calculada no es que sea mayor que 1, sino que le da una pendiente logarítmica de aproximadamente .5 con k, mientras que en el papel la pendiente es de aproximadamente 1.2, por lo tanto obtendrá la epsilon mal .

Mi algoritmo

Comencé copiando su código, revisándolo y volviéndolo a implementar de una manera equivalente. Mi reimplementación replicó sus resultados. Estoy bastante seguro de que implementó la versión discreta de la fórmula para E_b (k) correctamente. Sin embargo, un examen detallado del artículo sugiere que los autores usaron aproximaciones suaves en su código.

En la segunda página y columna, se establece la igualdad P (k | k '') = P (k, k'') / (k '') ^ (1-gamma). Esto es equivalente a reemplazar la probabilidad exacta P (k '') en el denominador de la primera integral con la aproximación suave de la ley de potencia (k'') ^ (- gamma) de la distribución de grados, y no es una igualdad.

El hecho de que los autores declaren esta aproximación como una igualdad sin calificación me sugiere que pueden haberla utilizado como tal en su código. Entonces, decidí usar su aproximación en el código, dando como resultado lo siguiente (donde obtuve gamma = 2.8 para cond-mat se explica a continuación).

def ebkss(g, b, gamma=2.8):
    edge_dict = defaultdict(lambda: defaultdict(int))
    degree_dict = defaultdict(int)
    edge_degree = [sorted(g.degree(e).values()) for e in g.edges()]
    for e in edge_degree:
        edge_dict[e[0]][e[-1]] +=1
    for i in g.degree().values():
        degree_dict[i] +=1
    edge_number = g.number_of_edges()
    node_number = g.number_of_nodes()
    ebks, ks = [], []
    for k1 in edge_dict:
        p1, p2 = 0, 0
        nk2k = np.sum(edge_dict[k1].values())
        pk1 = float(degree_dict[k1])/node_number
        k1pk1 = k1*pk1

        for k2 in edge_dict[k1]:
            if k2 >= b*k1:
                pk2k = float(edge_dict[k1][k2])/edge_number
                pk2 = float(degree_dict[k2])/node_number
                p1 += pk2k/(k2*k2**(-gamma))
        for k in degree_dict:
            if k>=b*k1:
                pk = float(degree_dict[k])/node_number
                p2 += pk
        if p2 > 0 and p1 > 0:
            ebks.append(p1/p2)
            ks.append(k1)
    return ebks, ks

Los resultados

Utilizando este código:

def get_logslope(x,y):
    A = np.empty((len(x), 2))
    A[:,0] = np.log(x)
    A[:,1] = 1
    res = la.lstsq(A, np.log(y))
    return res[0]

def show_eb(ca, b, gamma):
    #calculate ebk 
    ebk, k = ebkss(ca, b=b,gamma=gamma)
    print "Slope = ", get_logslope(np.array(k), np.array(ebk) )
    plt.plot(k,ebk,''r^'')
    plt.xlabel(r''$k$'', fontsize = 16)
    plt.ylabel(r''$E_b(k)$'', fontsize = 16)
    plt.xscale(''log'')
    plt.yscale(''log'')
    plt.show()
show_eb(ca, 3, 2.8)

Tengo esta salida:

Slope =  1.22136715547

La pendiente (hasta 1 dígito después del punto decimal, que es todo lo que se da en el documento) es correcta, y por lo tanto, épsilon ahora se puede calcular correctamente.

Sobre gamma

Obtuve el valor de gamma = 2.8 al agregar la pendiente de 1.2 al valor épsilon de 1.6 (esto se deduce de la ecuación 4 del documento). También realicé una comprobación rápida de la cordura con el módulo Powerlaw Python para determinar si esta gamma era adecuada.

import powerlaw
res = powerlaw.Fit(np.array(ca.degree().values())+1, xmin=10)
print res.alpha

Esta salida

2.84571139756

por lo tanto, 2.8 es correcto para el valor de gamma hasta el redondeo.

Editar con datos WWW

Probé mi método con el conjunto de datos WWW. Terminé obteniendo una pendiente que era cercana a la del papel, pero la escala aún está desactivada. Aquí está mi código:

def log_binning(x, y, bin_count=50):
    max_x = np.log10(max(x))
    max_y = np.log10(max(y))
    max_base = max([max_x,max_y])
    xx = [i for i in x if i>0]
    min_x = np.log10(np.min(xx))
    bins = np.logspace(min_x,max_base,num=bin_count)
    hist = np.histogram(x,bins)[0]
    nonzero_mask = np.logical_not(hist==0)       
    hist[hist==0] = 1
    bin_means_y = (np.histogram(x,bins,weights=y)[0] / hist)
    bin_means_x = (np.histogram(x,bins,weights=x)[0] / hist)
    return bin_means_x[nonzero_mask],bin_means_y[nonzero_mask]
def single_line_read(fname):    
    g = nx.Graph()
    with open(fname, "r") as f:
        for line in f:
          a = map(int,line.strip().split(" "))
          g.add_edge(a[0], a[1])
    return g

www = single_line_read("data/www.dat")
ebk, k = ebkss(www, 3, 2.6)
lk, lebk = log_binning(np.array(k,dtype=np.float64), np.array(ebk), bin_count=70)
#print lk, lebk
print "Slope", get_logslope(lk, lebk)
plt.plot(lk,lebk/www.number_of_edges(),''r^'')
plt.xlabel(r''$k$'', fontsize = 16)
plt.ylabel(r''$E_b(k)$'', fontsize = 16)
plt.xscale(''log'')
plt.yscale(''log'')
plt.show()

Pendiente 0.162453554297

La pendiente del papel original es 0.15. Obtuve el valor gamma de 2.6 mirando la Figura 3 en el documento (el gráfico gamma-épsilon).

En conclusión

No estoy seguro de por qué Eb (k) es mucho más pequeño que 1 en el gráfico del documento. Estoy bastante seguro de que se están realizando algunos ajustes de escala que no están explícitos en el documento. Sin embargo, pude recuperar el valor correcto de épsilon utilizando Eb (k). Siempre y cuando puedas calcular correctamente el épsilon, no me preocuparía demasiado.

Answer 3

Teniendo en cuenta el uso de la combinación de registros de los datos, se puede adoptar la siguiente función.

import numpy as np

def log_binning(x, y, bin_count=35):
    max_x = np.log10(max(x))
    max_y = np.log10(max(y))
    max_base = max([max_x,max_y])
    xx = [i for i in x if i>0]
    min_x = np.log10(np.min(xx))
    bins = np.logspace(min_x,max_base,num=bin_count)
    bin_means_y = (np.histogram(x,bins,weights=y)[0] / np.histogram(x,bins)[0])
    bin_means_x = (np.histogram(x,bins,weights=x)[0] / np.histogram(x,bins)[0])
    return bin_means_x,bin_means_y

Si desea agrupar linealmente los datos, use la siguiente función:

def LinearBinData(x, y, number): 
    data=sorted(zip(x,y))
    rs = np.linspace(min(x),max(x),number)
    rs = np.transpose(np.vstack((rs[:-1],rs[1:])))
    ndata = []
    within = []
    for start,end in rs:
        for i,j in data:
            if i>=start and i<end:
                within.append(j)
        ndata.append([(start+end)/2.0,np.mean(np.array(within))]  )
    nx,ny = np.array(ndata).T
    return nx,ny

Por lo general, para la relación de escalado, la combinación de registros sería una mejor opción.