2D-StruInforEntropy（Graph 二维结构信息熵）-- 记录

import numpy as np
import networkx as nx
import community
import matplotlib.pyplot as plt

# 无权重图 or all weights == 1.0
G = nx.Graph()
edges_list = [(0, 1), (0, 2), (0, 3), (0, 4), (1, 2), (1, 3), (2, 3),
              (4, 5), (4, 6), (4, 7), (4, 8), (5, 6), (5, 7), (5, 8), (6, 7), (6, 8), (7, 8), (7, 10),
              (9, 10), (9, 12), (9, 13), (10, 11), (10, 12), (11, 12), (11, 13), (12, 13)]
G.add_edges_from(edges_list)

G = nx.truncated_cube_graph()# eq
G = nx.truncated_tetrahedron_graph()# less
G = nx.tutte_graph()# better
G = nx.krackhardt_kite_graph()# less
G = nx.icosahedral_graph()# less
G = nx.house_x_graph()# better
G = nx.frucht_graph()# eq
G = nx.chvatal_graph()# # eq
G = nx.generators.barabasi_albert_graph(n=30, m=15)# less
G = nx.random_regular_graph(d=6, n=30)# less
G = nx.connected_watts_strogatz_graph(n=30, k=5, p=0.15)# better
# G = nx.grid_2d_graph(6, 6)# less
# G = nx.generators.barabasi_albert_graph(n=80, m=60)# less

nx.draw(G, node_size=500, with_labels=True, node_color=‘red‘)
plt.show()

nodes = [node for node in G.nodes()]
# G_du = G.degree()
# G_volume = sum(G_du[index_node] for index_node in G.nodes)
# G_pu_dic = {index_node: G_du[index_node]*1.0/(1.0*G_volume) for index_node in G.nodes}
# G_pu = [G_pu_dic[index_node] for index_node in G.nodes]

def pLog2p(p_i, eps = 1e-10):
    ind = np.where(p_i < eps)
    if len(ind) > 0:
        p_i[p_i < eps] = 1.0
        return p_i*np.log2(p_i)
    else:
        return p_i*np.log2(p_i)

def Log2p(p_i, eps = 1e-10):
    ind = np.where(p_i < eps)
    if len(ind) > 0:
        p_i[p_i < eps] = 1.0
        return np.log2(p_i)
    else:
        return np.log2(p_i)

def getGraphFractalDimension(G, plotFlg=False):
    from box_cover import GC
    X, Y = GC(G)
    LogXI = [np.log10(x) for x in X]
    LogYI = [np.log10(y) for y in Y]
    z1 = np.polyfit(LogXI, LogYI, 1)# 最小二乘多项式拟合
    fractal_dimension = z1[0]# The fractal dimension is the slope of the fitting line z1. 拟合线的斜率

    if plotFlg:
        import matplotlib.pyplot as plt
        p1 = np.poly1d(z1)  # 生成一个多项式
        yvals = p1(LogXI)
        plt.plot(LogXI, LogYI, ‘o‘, markersize=15)
        plt.plot(LogXI, yvals, ‘-‘, markersize=15, linewidth=3)
        plt.show()

    return fractal_dimension

def get_oneStruInforEntropy(G, weight=None):
    G_du = G.degree(weight=weight)
    G_volume = sum(G_du[index_node] for index_node in G.nodes)
    G_pu_dic = {index_node: G_du[index_node]*1.0/(1.0*G_volume) for index_node in G.nodes}
    G_pu = [G_pu_dic[index_node] for index_node in G.nodes]
    # Shonnon Entropy
    HP_G_Shonnon = sum(pLog2p(np.array(G_pu)))*(-1.0)
    return HP_G_Shonnon

# StruInforEntropy
def StruInforEntropy(G, partition):
    # nodes = G.nodes
    # n = G.number_of_nodes()
    sub_set = partition.copy()
    m = G.number_of_edges()
    degree = G.degree()
    Vij, gij, deg_ij = [], [], []
    for ind in range(len(sub_set)):
        sub_degree = 0
        dij = []
        for node in sub_set[ind]:
            sub_degree += degree[node]
            dij.append(degree[node])
        gj_c = nx.cut_size(G, sub_set[ind])
        Vij.append(sub_degree)
        gij.append(gj_c)
        deg_ij.append(np.array(dij))
    gij = np.array(gij, dtype=np.float)
    Vij = np.array(Vij, dtype=np.float)
    # print(‘Vi:‘, Vij)
    # print(‘gi:‘, gij)
    deg_ij = np.array(deg_ij)
    p_i = np.array([deg_ij[i] / Vij[i] for i in range(len(Vij))])
    pLogp = np.array([pLog2p(pi, eps=1e-10) for pi in p_i])
    sum_pLogp = np.array([np.sum(plp) for plp in pLogp])
    first = np.sum((-1.0) * Vij / (2.0 * m) * sum_pLogp)
    second = np.sum((-1.0) * gij / (2.0 * m) * Log2p(Vij / (2.0 * m)))
    HG = first + second

    return HG, Vij, gij, deg_ij

# # test: 2D-StruInforEntropy
# partition = [{0, 1, 2, 3}, {4, 5, 6, 7, 8}, {9, 10, 11, 12, 13}]
# HP_G, Vj, g_j, Gj_deg = StruInforEntropy(G, partition)

def delt_Xij(Xi, Xj, G):
    Xij = Xi + Xj
    sub_set = [Xi, Xj, list(set(Xij))]
    # sub_set = [[0, 1], [2, 3], [0, 1, 2, 3]]
    # n = G.number_of_nodes()
    m = G.number_of_edges()
    degree = G.degree()
    Vij, gij = [], []
    for ind in range(len(sub_set)):
        sub_degree = 0
        for node in sub_set[ind]:
            sub_degree += degree[node]
        gj_c = nx.cut_size(G, sub_set[ind])
        Vij.append(sub_degree)
        gij.append(gj_c)
    gij = np.array(gij)
    Vij = np.array(Vij)
    g_i, g_j, g_ij = gij[0], gij[1], gij[2]
    V_i, V_j, V_ij = Vij[0], Vij[1], Vij[2]
    log_Vij = Log2p(Vij, eps=1e-10)
    delt_G_Pij = 1.0 / 2.0 * m * ((V_i - g_i) * log_Vij[0] +
                                  (V_j - g_j) * log_Vij[1] -
                                  (V_ij - g_ij) * log_Vij[2] +
                                  (g_i + g_j - g_ij) * np.log2(2.0 * m))
    return delt_G_Pij

def doWhile(NodeA):
    NodeB = NodeA.copy()
    count = 0
    L = len(NodeA)
    for Xi in NodeB:
        for Xj in NodeB:
            if Xi != Xj:
                delt_ij = delt_Xij(Xi, Xj, G)# 函数
                if delt_ij <= 0:
                    count += 1
    if count - L*(L-1) != 0:
        return True
    return False

# Input 算法主体 -- 开始
print("Partition by min2DHG ..........")
global NodeA
nodes = list(G.nodes())
nodes.sort()
NodeA = [[node] for node in nodes]
print("Init-Input:", NodeA)
doWhileFlg = True
NodeA.reverse()
while doWhileFlg:
    Xi = NodeA.pop()
    Nj = NodeA.copy()
    delt_max = 0
    Xj_m = None
    for Xj in Nj:
        delt_ij = delt_Xij(Xi, Xj, G)# 函数
        if delt_ij > 0 and delt_ij > delt_max:
            Xj_m = Xj
            delt_max = delt_ij
    if Xj_m in Nj and Xj_m is not None:
        Nj.remove(Xj_m)
        Xij = Xi + Xj_m
        Nj.insert(0, Xij)
        # print(‘Xi:‘, Xi, ‘+ Xj:‘, Xj_m, ‘-->‘, Xij, ‘ delt_ij_HG:‘, delt_max)
    elif Xj_m is None:
        Nj.insert(0, Xi)
    NodeA = Nj
    doWhileFlg = doWhile(NodeA)
    # print(NodeA)
print(‘Output:‘, NodeA)
sub_set = NodeA.copy()
# Output: NodeA 算法主体 -- 结束

# m = G.number_of_edges()
# degree = G.degree()
# Vij, gij, deg_ij = [], [], []
# for ind in range(len(sub_set)):
#     sub_degree = 0
#     dij = []
#     for node in sub_set[ind]:
#         sub_degree += degree[node]
#         dij.append(degree[node])
#     gj_c = nx.cut_size(G, sub_set[ind])
#     Vij.append(sub_degree)
#     gij.append(gj_c)
#     deg_ij.append(np.array(dij))
# gij = np.array(gij, dtype=np.float)
# Vij = np.array(Vij, dtype=np.float)
# print(‘Vi:‘, Vij)
# print(‘gi:‘, gij)
# deg_ij = np.array(deg_ij, dtype=object)
# p_i = np.array([deg_ij[i]/Vij[i] for i in range(len(Vij))], dtype=object)
# pLogp = np.array([pLog2p(pi, eps=1e-10) for pi in p_i], dtype=object)
# sum_pLogp = np.array([np.sum(plp) for plp in pLogp])
# first = np.sum((-1.0) * Vij / (2.0 * m) * sum_pLogp)
# second = np.sum((-1.0) * gij / (2.0 * m) * Log2p(Vij / (2.0 * m)))
# HG = first + second
# print("Partition:", NodeA)
# print("StruInforEntropy:", HG)

# 2D-StruInforEntropy
HG, Vj, g_j, Gj_deg = StruInforEntropy(G, partition=sub_set)
print("Partition-Size(min2DHG):", len(sub_set))
print("Partition(min2DHG):", sub_set)
print("StruInforEntropy:", HG)

print("Partition by community ..........")
partition = community.best_partition(G)
size = float(len(set(partition.values())))
result = []
count = 0.
for com in set(partition.values()):
    count = count + 1.
    list_nodes = [nodes for nodes in partition.keys() if partition[nodes] == com]
    result.append(list_nodes)
print("Partition-Size(louvain):", len(result))
print("Partition(louvain):", result)
HG, Vj, g_j, Gj_deg = StruInforEntropy(G, partition=result)
print("StruInforEntropy:", HG)

print("Partition by pymetis ..........")
import pymetis
adjacency_list = [np.array([nei for nei in G.neighbors(n=node)]) for node in G.nodes()]
# adjacency_list = [np.array([4, 2, 1]),
#                   np.array([0, 2, 3]),
#                   np.array([4, 3, 1, 0]),
#                   np.array([1, 2, 5, 6]),
#                   np.array([0, 2, 5]),
#                   np.array([4, 3, 6]),
#                   np.array([5, 3])]
nparts = len(sub_set)
n_cuts, membership = pymetis.part_graph(nparts=nparts, adjacency=adjacency_list)
result = []
for npa in range(nparts):
    nodes_part = list(np.argwhere(np.array(membership) == npa).ravel())
    result.append(nodes_part)
print("Partition-Size(pymetis):", len(result))
print("Partition(pymetis):", result)
HG, Vj, g_j, Gj_deg = StruInforEntropy(G, partition=result)
print("StruInforEntropy:", HG)