import numpy as np import networkx as nx import infomap import community import matplotlib.pyplot as plt import matplotlib.colors as colors # 无权重图 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.karate_club_graph() 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) print(infomap.Infomap().version) def findCommunities(G): """ Partition network with the Infomap algorithm. Annotates nodes with 'community' id and return number of communities found. """ infomapX = infomap.Infomap("--two-level") # print("Building Infomap network from a NetworkX graph...") for e in G.edges(): infomapX.addLink(*e) print("Find communities with Infomap...") infomapX.run() print("Found {} modules with codelength: {}".format(infomapX.numTopModules(), infomapX.codelength)) communities = {} for node in infomapX.iterLeafNodes(): communities[node.physicalId] = node.moduleIndex() nx.set_node_attributes(G, values=communities, name='community') return infomapX.numTopModules(), communities def drawNetwork(G): # position map pos = nx.spring_layout(G) # community ids communities = [v for k, v in nx.get_node_attributes(G, 'community').items()] numCommunities = max(communities) + 1 # color map from http://colorbrewer2.org/ cmapLight = colors.ListedColormap(['#a6cee3', '#b2df8a', '#fb9a99', '#fdbf6f', '#cab2d6'], 'indexed', numCommunities) cmapDark = colors.ListedColormap(['#1f78b4', '#33a02c', '#e31a1c', '#ff7f00', '#6a3d9a'], 'indexed', numCommunities) # Draw edges nx.draw_networkx_edges(G, pos) # Draw nodes nodeCollection = nx.draw_networkx_nodes(G, pos = pos, node_color = communities, cmap = cmapLight) # Set node border color to the darker shade darkColors = [cmapDark(v) for v in communities] nodeCollection.set_edgecolor(darkColors) # Draw node labels for n in G.nodes(): plt.annotate(n, xy = pos[n], textcoords = 'offset points', horizontalalignment = 'center', verticalalignment = 'center', xytext = [0, 0], color = cmapDark(communities[n])) plt.axis('off') # plt.savefig("karate.png") plt.show() numModules, communities = findCommunities(G) drawNetwork(G) list_values = [val for val in communities.values()] list_nodes = [node for node in communities.keys()] result = [] for ndpa in range(numModules): nodes_part_index = np.argwhere(np.array(list_values) == ndpa).ravel() nodes_part = list(np.array(list_nodes)[nodes_part_index]) result.append(nodes_part) print("Partition-Size(Infomap):", len(result)) print("Partition(Infomap):", result) HG, Vj, g_j, Gj_deg = StruInforEntropy(G, partition=result) print("StruInforEntropy:", HG)
2D-StruInforEntropy（Graph 二维结构信息熵）-- 记录
个人学习记录