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

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

思考

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

文中

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

思考

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


import numpy as np
import matplotlib.pyplot as pltclass SIR(object):def __init__(self, beta, gama, s0, i0, r0):"""初始化模型参数:param beta: 感染系数:param gama: 恢复系数"""self.beta = betaself.gama = gamaself.initalization(s0, i0, r0)def initalization(self, s0, i0, r0):"""初始化易感者、感染者和移除者在总体人数N中的比例:param s0::param i0::param r0::return:"""self.s0 = s0self.i0 = i0self.r0 = r0def sir_model(self, y, beta, gama):s, i, r = yds_dt = - beta * s * idi_dt = beta * s * i - gama * idr_dt = gama * ireturn np.array([ds_dt, di_dt, dr_dt])def calc(self):# 自定义计算拟合函数，不使用odeint函数，以方便后期的图上扩展resu = [] # 计算结果x_ = []  # x坐标轴，存储循环计算次数x = 0while abs(self.i0) > pow(10, -6):x_.append(x)x += 1ds_dt, di_dt, dr_dt = self.sir_model([self.s0, self.i0, self.r0], self.beta, self.gama)self.s0 += ds_dtself.i0 += di_dtself.r0 += dr_dtresu.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 pdclass MCIM(object):def __init__(self, G, k):""":param G:  networkx:param k:  种子节点大小"""self.G = G  # 图self.k = kself.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 = 0for 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 * resultdef calcNodeIDS_E2_dv2(self, node):resu = 0for u in list(self.G.neighbors(node)):resu += self.G.degree(u)return resudef findMaxCalcNodeIDS_E2_dv2(self):max = 0for i in range(len(self.G)):val = self.calcNodeIDS_E2_dv2(i)if val > max:max = valreturn maxdef calcNodeIDS_Coefficient(self, v):return self.calcNodeIDS_E2_dv2(v) / self.findMaxCalcNodeIDS_E2_dv2()def calcNodeIDS_E2(self, v):result = 0for 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 * resultdef calcNodeDS(self, node):# calculate node v's degreereturn 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 = 0for u in list(self.G.neighbors(node)):if u in self.SS:resu += 1return resudef calcNodeIDO(self, v):Nv = self.G.neighbors(v)resu = 0for seed in self.SS:Nu = self.G.neighbors(seed)com = [val for val in Nv if val in Nu]temp_resu = 0for node in com:if node in self.SS:temp_resu += 1resu += temp_resureturn resu# 生成矩阵Adef 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中的数据类型dataframedata1 = 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] += 1for 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)# endreturn 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 pltclass SIR(object):def __init__(self, beta, gama, s0, i0, r0, G, seed_node_list):"""初始化模型参数:param beta: 感染系数:param gama: 恢复系数"""self.beta = betaself.gama = gamaself.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 = s0self.i0 = i0self.r0 = r0self.G = Gself.seed_node_list = seed_node_listdef sir_model(self, y, beta, gama):s, i, r = yds_dt = - beta * s * idi_dt = beta * s * i - gama * idr_dt = gama * ireturn np.array([ds_dt, di_dt, dr_dt])def calc(self):# 自定义resu = [] # 计算结果x_ = []  # x坐标轴，存储循环计算次数x = 0influenced_node = []G = self.Gi_state_node = self.seed_node_list[:] # 列表拷贝r_state_node = []count = 0count_ = []resu = []# 假设每个节点之间的传染行为是相互独立的# 对于每个节点来说，传播的初始s/i/r是累计的s0, i0, r0 = self.s0, self.i0, self.r0while 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 += 1count_.append(count)resus = self.sir_model([s0, i0, r0], beta, self.gama)# 更新s0/i0/r0s0 += resus[0]i0 += resus[1]r0 += resus[2]resu.append([s0, i0, r0])# 该节点状态的改变m = 0for v in G.neighbors(u):if v in i_state_node:m += 1pro = 1 - pow(1 - beta, m)if random.random() < pro:# 节点状态变为Ii_state_node.append(u)influenced_node.append(u)# 当前节点node也有一定的几率变为Rif random.random() < self.gama:# 节点状态变为Rr_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 MCIMG = 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 pltclass draw(object):def __init__(self, G, seed_set):self.G = Gself.seeds = seed_setdef initialization(self):others = []labels = {}pos = nx.spring_layout(self.G)  # positions for all nodesfor i in range(len(self.G)):labels[i] = str(i)if i not in self.seeds:others.append(i)return others, labels, posdef 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# nodesnx.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)# edgesnx.draw_networkx_edges(self.G, pos, width=1.0, alpha=0.5)# labelsnx.draw_networkx_labels(self.G, pos, labels, font_size=14, font_color='white')# showplt.axis("off") # 关闭坐标轴plt.show()# 程序入口
if __name__ == '__main__':import MCIM, SIRk = 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