SIR模型的应用（2） - Influence maximization in social networks based on TOPSIS(3)

上篇中我们说道了，需要注意作者是如何将SIR应用到种子节点的传播过程之中的。那么我们就来看看作者是如何做的。

思考

不妨先思考下： ① 初始时刻，处于I状态的是种子节点，其余节点是易感节点； ② 在进行传播时，传播的途径是节点之间的连接边；也就是在SIR模型进行传播的时候，需要满足两个条件： 1）该节点到目标节点之间有直接连接边； 2）待传播节点是易感节点，即S；

文中

均匀模型，图中所有边的传播概率均是β；基于度的模型，节点v到节点u的传播概率可以计算为$1/d_u$

思考

由于，我们在SIR模型中需要对图中的节点进行计算，故而就不是简单的用odeint函数进行求解拟合。所以，这里需要自己写一个函数来模拟计算odeint函数的效果，才能在后面的步骤中加入对单个节点的考虑。

import numpy as np import matplotlib.pyplot as plt class SIR(object): def __init__(self, beta, gama, s0, i0, r0): """ 初始化模型参数 :param beta: 感染系数 :param gama: 恢复系数 """ self.beta = beta self.gama = gama self.initalization(s0, i0, r0) def initalization(self, s0, i0, r0): """ 初始化易感者、感染者和移除者在总体人数N中的比例 :param s0: :param i0: :param r0: :return: """ self.s0 = s0 self.i0 = i0 self.r0 = r0 def sir_model(self, y, beta, gama): s, i, r = y ds_dt = - beta * s * i di_dt = beta * s * i - gama * i dr_dt = gama * i return np.array([ds_dt, di_dt, dr_dt]) def calc(self): # 自定义计算拟合函数，不使用odeint函数，以方便后期的图上扩展 resu = [] # 计算结果 x_ = [] # x坐标轴，存储循环计算次数 x = 0 while abs(self.i0) > pow(10, -6): x_.append(x) x += 1 ds_dt, di_dt, dr_dt = self.sir_model([self.s0, self.i0, self.r0], self.beta, self.gama) self.s0 += ds_dt self.i0 += di_dt self.r0 += dr_dt resu.append([self.s0, self.i0, self.r0]) return np.array(x_), np.array(resu) def plot(self): time, res = self.calc() plt.figure() plt.plot(time, res[:, 0], label="S(t)") plt.plot(time, res[:, 1], label="I(t)") plt.plot(time, res[:, 2], label="R(t)") plt.legend() plt.grid() plt.xlabel("计算次数") plt.ylabel("proportions") plt.title("SIR model simulation") plt.show()

不难发现，上面的calc函数中，计算也可以拟合计算SIR模型的计算过程。那么，我们接着就考虑如何将这个模型应用到图上。

在基于均匀模型中，图中所有边的传播概率均是β，文中分别探究了β取值对传播的影响，如下图：

文中的实验中，β取值为0.05，由于文中同样对γ进行了实验，如下图： 非常有意思，就是当γ取值比较大的时候，SIR模型的图像，就不是那么规则了，如下图γ取值2.2：

由于我们判断的终止条件定义为处于感染状态的节点全部消失，那么在图中，当没有与之相连的S状态的节点的时候，模型中的S到R的计算应该继续。可以简单的将此时的β定义为0。

实验

不妨采用基于度的模型进行传播，将该值作为β值，而γ的取值，就取0.05。进行编写在图上的实验。 初始选定6个节点作为种子节点：

使用前面文中中实现的MCIM算法来进行种子节点集的选择，编写基于图结构传播的的SIR算法，来模拟病毒的传播过程，传染过的节点如下： 不妨打印下种子节点、处于R状态的节点、被影响节点： 完整程序分为种子节点选择MCIM.py、SIR.py、DrawGraph.py三个文件，分别如下：

1. MCIM.py

# Influence maximization in social networks based on TOPSIS import networkx as nx import math import numpy as np import pandas as pd class MCIM(object): def __init__(self, G, k): """ :param G: networkx :param k: 种子节点大小 """ self.G = G # 图 self.k = k self.SS = [] # 种子节点集合 self.SS.append(self.getMaxDegreeNode()) # 初始加入最大度节点 self.weight = [0.3, 0.2, 0.3, 0.2] # topsis判别计算权重 # 得到图中度最大的节点 def getMaxDegreeNode(self): dic = self.G.degree() return sorted(dict(dic).items(), key=lambda x: x[1], reverse=True)[:1][0][0] def calcNodeIDS_E1(self, v): result = 0 for u in list(self.G.neighbors(v)): result += self.G.degree(u) / self.G.degree(v) * math.log(self.G.degree(u) / self.G.degree(v)) return -1 * result def calcNodeIDS_E2_dv2(self, node): resu = 0 for u in list(self.G.neighbors(node)): resu += self.G.degree(u) return resu def findMaxCalcNodeIDS_E2_dv2(self): max = 0 for i in range(len(self.G)): val = self.calcNodeIDS_E2_dv2(i) if val > max: max = val return max def calcNodeIDS_Coefficient(self, v): return self.calcNodeIDS_E2_dv2(v) / self.findMaxCalcNodeIDS_E2_dv2() def calcNodeIDS_E2(self, v): result = 0 for u in list(self.G.neighbors(v)): result += self.G.degree(u) / self.calcNodeIDS_E2_dv2(v) * math.log( self.G.degree(u) / self.calcNodeIDS_E2_dv2(v)) return -1 * result def calcNodeDS(self, node): # calculate node v's degree return self.G.degree(node) def calcNodeIDS(self, node): return self.calcNodeIDS_E1(node) + self.calcNodeIDS_Coefficient(node) * self.calcNodeIDS_E2(node) def calcNodeDO(self, node): resu = 0 for u in list(self.G.neighbors(node)): if u in self.SS: resu += 1 return resu def calcNodeIDO(self, v): Nv = self.G.neighbors(v) resu = 0 for seed in self.SS: Nu = self.G.neighbors(seed) com = [val for val in Nv if val in Nu] temp_resu = 0 for node in com: if node in self.SS: temp_resu += 1 resu += temp_resu return resu # 生成矩阵A def generateMatrixA(self): A = [] for i in range(len(self.G)): # 依次计算每个节点的四个特征 A.append([self.calcNodeIDS(i), self.calcNodeIDS(i), -1 * self.calcNodeDO(i), -1 * self.calcNodeIDO(i)]) return np.array(A) # 决策函数 def topsis(self, matrixa, weight, n): # 二维数组的矩阵转换成pandas中的数据类型dataframe data1 = pd.DataFrame(matrixa, index=[i for i in range(n)], columns=['DS', 'IDS', 'DO', 'IDO']) # 归一化 data = data1 / np.sqrt((data1 ** 2).sum()) # 最优最劣方案 # Z是正理想解和负理想解矩阵 Z = pd.DataFrame([data.min(), data.max()], index=['负理想解', '正理想解']) # 距离 # weight = entropyWeight(data) if weight is None else np.array(weight) Result = data1.copy() Result['正理想解'] = np.sqrt(((data - Z.loc['正理想解']) ** 2 * weight).sum(axis=1)) Result['负理想解'] = np.sqrt(((data - Z.loc['负理想解']) ** 2 * weight).sum(axis=1)) # 综合得分 Result['综合得分'] = Result['负理想解'] / (Result['负理想解'] + Result['正理想解']) Result['排序'] = Result.rank(ascending=False, method='first')['综合得分'] return Result[(Result['排序'] == 1.0)].index.tolist() def calc(self): MatrixA = self.generateMatrixA() removed_node_list = [] for i in range(self.k - 1): # remove u from matrix ， 即：删除所在行 np.delete(MatrixA, i, 0) node_i = self.SS[len(self.SS) - 1] removed_node_list.append(node_i) for node in list(self.G.neighbors(node_i)): if node not in removed_node_list: MatrixA[node][2] += 1 for i_node in list(self.G.neighbors(node)): if i_node not in removed_node_list: com = [val for val in list(self.G.neighbors(node_i)) if val in list(self.G.neighbors(node))] MatrixA[i_node][3] += len(com) # find the optimal solution by topsis(A) u = self.topsis(MatrixA, self.weight, len(MatrixA))[0] self.SS.append(u) # end return self.SS # 调用案例 if __name__ == '__main__': G = nx.karate_club_graph() print(MCIM(G, 5).calc())

2. SIR.py

import numpy as np import networkx as nx import random import matplotlib.pyplot as plt class SIR(object): def __init__(self, beta, gama, s0, i0, r0, G, seed_node_list): """ 初始化模型参数 :param beta: 感染系数 :param gama: 恢复系数 """ self.beta = beta self.gama = gama self.initalization(s0, i0, r0, G, seed_node_list) def initalization(self, s0, i0, r0, G, seed_node_list): """ 初始化易感者、感染者和移除者在总体人数N中的比例 :param s0: :param i0: :param r0: :param G: networkx 生成的图数据G :param G: list 种子节点列表 :return: """ self.s0 = s0 self.i0 = i0 self.r0 = r0 self.G = G self.seed_node_list = seed_node_list def sir_model(self, y, beta, gama): s, i, r = y ds_dt = - beta * s * i di_dt = beta * s * i - gama * i dr_dt = gama * i return np.array([ds_dt, di_dt, dr_dt]) def calc(self): # 自定义 resu = [] # 计算结果 x_ = [] # x坐标轴，存储循环计算次数 x = 0 influenced_node = [] G = self.G i_state_node = self.seed_node_list[:] # 列表拷贝 r_state_node = [] count = 0 count_ = [] resu = [] # 假设每个节点之间的传染行为是相互独立的 # 对于每个节点来说，传播的初始s/i/r是累计的 s0, i0, r0 = self.s0, self.i0, self.r0 while len(i_state_node) > 0: # 取出每个节点，进行SIR传播 node = i_state_node.pop() # 找邻居节点 node_neighbors = G.neighbors(node) # 依次计算每个邻居节点的度 for u in node_neighbors: # 判断该节点的状态 S/I/R # 邻居节点处于以感染状态，先不计算它 # 邻居节点处于以恢复状态，先不计算它 if not (u in i_state_node or u in r_state_node): # 邻居节点处于S状态 # 计算感染系数 beta = 1 / G.degree(u) # 进行SIR模型的计算 if abs(i0) > pow(10, -6) > 0: count += 1 count_.append(count) resus = self.sir_model([s0, i0, r0], beta, self.gama) # 更新s0/i0/r0 s0 += resus[0] i0 += resus[1] r0 += resus[2] resu.append([s0, i0, r0]) # 该节点状态的改变 m = 0 for v in G.neighbors(u): if v in i_state_node: m += 1 pro = 1 - pow(1 - beta, m) if random.random() < pro: # 节点状态变为I i_state_node.append(u) influenced_node.append(u) # 当前节点node也有一定的几率变为R if random.random() < self.gama: # 节点状态变为R r_state_node.append(node) else: # 节点保持I状态 i_state_node.append(node) # print("种子节点集合：", self.seed_node_list, ",len=", len(self.seed_node_list)) # print("影响的节点集合为：", [x for x in r_state_node if x not in self.seed_node_list], ",len=", len([x for x in r_state_node if x not in self.seed_node_list])) return r_state_node,influenced_node # 调用案例 if __name__ == '__main__': import MCIM G = nx.karate_club_graph() seed = MCIM.MCIM(G, 5).calc() SIR(None, 0.05, 1 - (5 / len(G)), 5 / len(G), 0, G, seed).calc()

3. DrawGraph.py

import networkx as nx import matplotlib.pyplot as plt class draw(object): def __init__(self, G, seed_set): self.G = G self.seeds = seed_set def initialization(self): others = [] labels = {} pos = nx.spring_layout(self.G) # positions for all nodes for i in range(len(self.G)): labels[i] = str(i) if i not in self.seeds: others.append(i) return others, labels, pos def draw(self): options = {"node_size": 500, "alpha": 0.8} others, labels, pos = self.initialization() # https://blog.csdn.net/zhou0107/article/details/9146299?utm_medium=distribute.pc_relevant.none-task-blog-BlogCommendFromMachineLearnPai2-3.add_param_isCf&depth_1-utm_source=distribute.pc_relevant.none-task-blog-BlogCommendFromMachineLearnPai2-3.add_param_isCf # nodes nx.draw_networkx_nodes(self.G, pos, nodelist=self.seeds, label=self.seeds, node_color="OrangeRed", **options) nx.draw_networkx_nodes(self.G, pos, nodelist=others, label=others, node_color="Silver", **options) # edges nx.draw_networkx_edges(self.G, pos, width=1.0, alpha=0.5) # labels nx.draw_networkx_labels(self.G, pos, labels, font_size=14, font_color='white') # show plt.axis("off") # 关闭坐标轴 plt.show() # 程序入口 if __name__ == '__main__': import MCIM, SIR k = 6 # 种子节点个数 G = nx.karate_club_graph() seed = MCIM.MCIM(G, k).calc() print("Seed_Nodes") print(seed) draw(G, seed).draw() r_state_node, influenced_node = SIR.SIR(None, 0.05, 1 - (k / len(G)), k / len(G), 0, G, seed).calc() print("R_State_Nodes") print(r_state_node) print("Influenced_Node") print(influenced_node) draw(G, seed + influenced_node).draw()

附：颜色名列表

不妨再改造改造程序，绘制一下再这个图上的SIR变化曲线如： 可以看出，和上篇《SIR模型的应用 - Influence maximization in social networks based on TOPSIS(3)》的SIR原始的图形基本是拟合的。

代码地址：here