Python大作业，题目如下，这里传下我的屎山代码以供参考。
目录
题目
代码
import sys
import numpy as np
import matplotlib.pyplot as plt
import networkx as nx
import random
from sklearn.metrics import f1_score, roc_auc_score
import warnings
import time
import multiprocessing as mp
np.set_printoptions(threshold=sys.maxsize) # 显示全部矩阵 调试更方便
warnings.filterwarnings('ignore') # 忽略futurewarning
plt.rcParams['font.sans-serif'] = ['SimHei'] # 中文字体 使得画图可以显示中文
def initialize_graph(N, p, num):
"""初始化num个节点为感染节点
:param N:网络规模
:param p:连边概率
:param num:初始感染数量
:return:(随机ER网络,初始感染节点列表)
"""
graph = nx.random_graphs.erdos_renyi_graph(N, p) # 生成包含N个节点、连边概率p的ER随机图
while not nx.is_connected(graph): # 重新生成直到是连通图 否则后面求逆会出错
graph = nx.random_graphs.erdos_renyi_graph(N, p)
# pos = nx.spring_layout(graph) # 节点分布形式 这里选取节点随机分布
# nx.draw_networkx(graph, pos, node_color="green", edge_color="red", with_labels=False, node_size=5)
# plt.show()
initial_infected = random.sample(list(graph.nodes()), num) # 随机选择num个节点作为初始感染节点
for node in graph.nodes(): # 遍历节点
if node in initial_infected: # 如果是初始感染节点 就把节点的status属性设为I 否则设为S
graph.nodes[node]["status"] = "I" # 初始感染节点状态为感染
else:
graph.nodes[node]["status"] = "S" # 其他节点初始状态为健康
return graph, sorted(initial_infected)
def SIR_model(graph, beta, miu, t_obs):
"""SIR模型 先感染再康复
:param graph:图
:param beta:感染概率
:param miu:康复概率
:param t_obs:观测时刻
:return:经过t_obs天后的图
"""
for _ in range(t_obs):
old_graph = graph.copy() # 备份
for node in old_graph.nodes():
if old_graph.nodes[node]["status"] == "S": # 如果是健康节点
neighbors = list(old_graph.neighbors(node)) # 找该健康节点的所有邻接节点
cnt = sum(1 for neighbor in neighbors if old_graph.nodes[neighbor]["status"] == "I") # 有cnt个感染节点邻居
if random.random() < 1 - (1 - beta) ** cnt: # 该节点被感染的概率为 1−(1−𝛽)^𝑚
graph.nodes[node]["status"] = "I" # 感染
for node in old_graph.nodes():
if old_graph.nodes[node]["status"] == "I": # 如果是感染节点
if random.random() < miu:
graph.nodes[node]["status"] = "R" # 康复
return graph
def cal_S(graph):
"""计算相似度矩阵S
:param graph:图
:return:
"""
adj_matrix = nx.adjacency_matrix(graph).toarray() # graph转邻接矩阵
row_sums = np.sum(adj_matrix, axis=1) # 求行和
D = np.diag(row_sums) # 转成对角矩阵
D_new = np.linalg.inv(np.sqrt(D)) # 求-0.5次幂 如果测试规模太小导致生成的图含有孤立结点 这里就会报错 因为含有孤立结点的图的邻接矩阵是个奇异矩阵 不能求逆
S = np.dot(np.dot(D_new, adj_matrix), D_new) # D^-0.5 * W * D^-0.5
return S
def label_rev(graph, S, alpha, max_iterations):
"""标签反向传播
:param graph:图
:param S:相似度矩阵
:param alpha:传播系数
:param max_iterations:最大迭代次数
:return:收敛或经过最大迭代次数的标签
"""
labels = [1 if graph.nodes[node]["status"] == "I" else -1 for node in graph.nodes()] # 感染状态标签设置为1 易感和康复设为-1
Y = np.array(labels).astype(float) # 初始标签向量Y 转成浮点型的array 不然是整型的
G = Y.copy() # 一开始为初始标签
for iter in range(max_iterations): # 迭代次数
old_G = G.copy() # 拷贝上一次迭代的标签 拿来判断收敛的
for i in range(len(Y)):
neighbors = list(graph.neighbors(i)) # 找所有邻接节点
# print("第{}次迭代 节点{}的邻居是{}".format(iter, i, neighbors))
label_sum = sum([S[i][j] * G[j] for j in neighbors]) # 遍历邻接节点 根据相似度矩阵加权求和
G[i] = alpha * label_sum + (1 - alpha) * Y[i] # 更新标签
if np.allclose(old_G, G, rtol=1e-16, atol=1e-16): # 收敛 old_G和G非常相似的话就判定收敛 停止迭代 相对容差和绝对误差这里设的很小 差不多是双精度极限了
# print("old_G={}\nG={}".format(old_G,G))
print("在第{}次迭代收敛".format(iter))
break
else:
print("到达最大迭代次数")
return G
def find_source(graph, G):
"""溯源找初始传播者 若节点标签大于所有邻接节点标签和 则认为该节点是初始传播者
:param graph:图
:param G:标签
:return:初始传播者列表
"""
predict_infected = []
for node in graph.nodes():
adj_matrix = nx.adjacency_matrix(graph).toarray() # graph转邻接矩阵
if G[node] > max(adj_matrix[node] * G): # 如果该节点大于所有邻居节点的标签 即认为是初始传播者
predict_infected.append(node)
return predict_infected
# return [node for node in graph.nodes() if G[node] >max(nx.adjacency_matrix(graph).toarray()[node] * G) ] # 精简版 虽然推导式被优化过 但可读性垃圾
def accuracy(init_infected, predict_infected, N):
# 调库算准确率
true_labels = [1 if node in init_infected else 0 for node in range(N)] # 真实值：初始感染者为1，其他为0
predicted_labels = [1 if node in predict_infected else 0 for node in range(N)] # 预测值：预测感染者为1，其他为0
f1 = f1_score(true_labels, predicted_labels) # 调库算f1
auc = roc_auc_score(true_labels, predicted_labels) # 调库算auc
return f1, auc
def single_simulate(N, p, alpha, beta, miu, t_obs, max_iterations, num_initial_infected):
# 单次模拟 返回单个f1和auc的值
graph, init_infected = initialize_graph(N, p, num_initial_infected) # 初始化网络
result_graph = SIR_model(graph, beta, miu, t_obs) # 在t_bos天观测节点情况
S = cal_S(result_graph) # 计算相似度矩阵
G = label_rev(graph, S, alpha, max_iterations) # 标签迭代
# print("最终标签:{}".format(G))
predict_infected = find_source(graph, G) # 溯源推测初始传播者
print("真实初始传播者:{}".format(init_infected))
print("推测初始传播者:{}".format(predict_infected))
f1, auc = accuracy(init_infected, predict_infected, N) # 计算f1和auc
print("f1-score:{:.2f}".format(f1))
print("auc:{:.2f}".format(auc))
return f1, auc
def plot(x, y0, y1, xlabel):
# 画双Y轴折线图 y0是f1 y1是auc
fig, ax1 = plt.subplots(figsize=(12.8, 7.2))
ax2 = ax1.twinx()
line1 = ax1.plot(x, y0, lw=3, ls='-', c='b', label="f1-score")
line2 = ax2.plot(x, y1, lw=3, ls='-', c='r', label="auc") # 设置横坐标标签
plt.xticks(x) # 设置横坐标刻度
plt.title("{}-f1与auc准确率关系图".format(xlabel)) # 设置图表标题
ax1.set_xlabel(xlabel) # 设置横坐标标签
ax1.set_ylabel("f1准确率", color='b') # 设置y0纵坐标标签
ax1.tick_params(axis='y', labelcolor='b')
ax2.set_ylabel("auc准确率", color='r') # 设置y1纵坐标标签
ax2.tick_params(axis='y', labelcolor='r')
lines = line1 + line2
plt.legend(lines, [l.get_label() for l in lines], loc='best') # 图例
fig.tight_layout()
fig.savefig("{}-f1,auc准确率关系图".format(xlabel)) # 保存
def param_optimize(N, p, alpha, beta, miu, t_obs, max_iterations, num_initial_infected, param, paramrange, maxtimes):
# 改变参数 优化 每个参数都模拟maxtimes次 计算得到的f1和auc取均值
f1_list = []
auc_list = []
if param == "p":
for p in paramrange:
f1 = auc = 0 # 该参数的f1和auc总和
for _ in range(maxtimes):
single_f1, single_auc = single_simulate(N, p, alpha, beta, miu, t_obs, max_iterations,
num_initial_infected)
f1 += single_f1
auc += single_auc
f1_list.append(f1 / maxtimes) # 取均值
auc_list.append(auc / maxtimes) # 取均值
return paramrange, f1_list, auc_list, "连边概率p"
elif param == "beta":
for beta in paramrange:
f1 = auc = 0 # 该参数的f1和auc总和
for _ in range(maxtimes):
single_f1, single_auc = single_simulate(N, p, alpha, beta, miu, t_obs, max_iterations,
num_initial_infected)
f1 += single_f1
auc += single_auc
f1_list.append(f1 / maxtimes) # 取均值
auc_list.append(auc / maxtimes) # 取均值
return paramrange, f1_list, auc_list, "感染率beta"
elif param == "alpha":
for alpha in paramrange:
f1 = auc = 0 # 该参数的f1和auc总和
for _ in range(maxtimes):
single_f1, single_auc = single_simulate(N, p, alpha, beta, miu, t_obs, max_iterations,
num_initial_infected)
f1 += single_f1
auc += single_auc
f1_list.append(f1 / maxtimes) # 取均值
auc_list.append(auc / maxtimes) # 取均值
return paramrange, f1_list, auc_list, "传播概率alpha"
elif param == "t_obs":
for t_obs in paramrange:
f1 = auc = 0 # 该参数的f1和auc总和
for _ in range(maxtimes):
single_f1, single_auc = single_simulate(N, p, alpha, beta, miu, t_obs, max_iterations,
num_initial_infected)
f1 += single_f1
auc += single_auc
f1_list.append(f1 / maxtimes) # 取均值
auc_list.append(auc / maxtimes) # 取均值
return paramrange, f1_list, auc_list, "观测时刻t_obs"
if __name__ == "__main__":
N = 100 # 规模
p = 0.08 # 连边概率
alpha = 0.1 # 传播概率
beta = 0.15 # 感染概率
miu = 0.1 # 康复概率
t_obs = 4 # 观测时刻
max_iterations = 100 # 标签迭代次数
num_initial_infected = 5 # 初始感染人数
# print(single_simulate(N, p, alpha, beta, miu, t_obs, max_iterations, num_initial_infected))
# 取10次均值
# start_t = time.perf_counter() # 计时
# plot(*param_optimize(N, p, alpha, beta, miu, t_obs, max_iterations, num_initial_infected, param="p",
# paramrange=np.arange(0.06, 0.32, 0.02), maxtimes=10))
# plot(*param_optimize(N, p, alpha, beta, miu, t_obs, max_iterations, num_initial_infected, param="beta",
# paramrange=np.arange(0.05, 1.05, 0.05), maxtimes=10))
# plot(*param_optimize(N, p, alpha, beta, miu, t_obs, max_iterations, num_initial_infected, param="alpha",
# paramrange=np.arange(0.1, 1.1, 0.1), maxtimes=10))
# plot(*param_optimize(N, p, alpha, beta, miu, t_obs, max_iterations, num_initial_infected, param="t_obs",
# paramrange=np.arange(2, 9, 1), maxtimes=10))
# end_t = time.perf_counter()
# print("消耗: " + "{:.2f}".format(end_t-start_t) + " 秒")
# 下面是多进程计算 自己的电脑用上面的单核计算72秒 下面并行计算36秒
start_t = time.perf_counter() # 计时
num_cores = int(mp.cpu_count())
print("本地计算机有: " + str(num_cores) + " 核心")
if num_cores <= 4:
pool = mp.Pool(num_cores)
else:
pool = mp.Pool(4) # 最多用四个核心
param_dict = {
'task1': (N, p, alpha, beta, miu, t_obs, max_iterations, num_initial_infected, "p",
np.arange(0.06, 0.32, 0.02),10),
'task2': (N, p, alpha, beta, miu, t_obs, max_iterations, num_initial_infected, "beta",
np.arange(0.05, 1.05, 0.05), 10),
'task3': (N, p, alpha, beta, miu, t_obs, max_iterations, num_initial_infected, "alpha",
np.arange(0.1, 1.1, 0.1),10),
'task4': (N, p, alpha, beta, miu, t_obs, max_iterations, num_initial_infected, "t_obs",
np.arange(2, 9, 1),10)}
results = [pool.apply_async(param_optimize, args=args) for args in param_dict.values()]
results = [p.get() for p in results]
# print(results)
for result in results:
plot(*result) # 画图
end_t = time.perf_counter()
print("多进程计算 消耗: " + "{:.2f}".format(end_t - start_t) + " 秒")
步骤
（偷懒直接复制报告里的内容了）
1. 生成ER随机网络：通过initialize_graph函数，生成了一个包含N个节点、连边概率为p的ER随机图，一开始用的列表来做整体的存储结构，取上半矩阵转置求和，后来发现有个networkx库，直接调用nx.random_graphs.erdos_renyi_graph(N, p)来生成ER随机网络了，值得注意的是这里要重复生成ER图直到是连通图，不然后面求相似度矩阵中求逆会报差错（奇异矩阵不能求逆）。初始时，随机选择了num个节点作为感染节点（后来取的是五个初始感染节点）。并把该节点的status属性改为I，代表感染节点，其余节点status属性改为S，代表健康节点。返回ER图和初始感染者列表。
2. SIR传播模型：利用SIR_model函数，实现了离散的SIR传播动力学模型。该模型在t_obs天观测节点状态，其中感染概率为beta，康复概率为miu。值得注意的是这里是先感染再康复。遍历t_obs，每次遍历时先拷贝原来的ER图（之前没备份导致出来的预测率很低，找了一万年bug）。模拟感染过程，遍历节点，如果是健康节点，计算邻居节点的感染数量，若随机数小于1−(1−𝛽)^𝑚则判定感染，status改为I。模拟康复过程，遍历节点，如果是感染节点，若随机数小于μ则判定康复，status改为R。值得注意的是，这里的康复节点假设已经有抗体，不会再被感染。
3. 相似度矩阵计算： cal_S函数,通过nx.adjacency_matrix将ER图转为邻接矩阵，计算每行的和，再将行和转对角矩阵，再求-0.5次幂（图含孤立结点就会报错）得到D^-0.5，最后返回D^-0.5 * W * D^-0.5，即为相似度矩阵S。
4. 标签反向传播：label_rev函数，进行标签反向传播迭代，推测感染源。Y是初始标签向量（这里采取了float型，精度越高越好），G一开始等于初始标签向量，在每次迭代中，先拷贝G为old_G，再遍历所有节点，计算该节点的所有邻居节点与相似度矩阵的乘积的和，label_sum = sum([S[i][j] * G[j] for j in neighbors])即公式
再更新标签G，G[i] = alpha * label_sum + (1 - alpha) * Y[i]。每次迭代后都要利用np.allclose(old_G, G, rtol=1e-16, atol=1e-16)判定old_G和G是不是相似，若绝对容差小于1e-16或者相对容差小于1e-16则判定收敛跳出大循环（这里给的容差很小，不然迭代个几次就收敛了）。或者到达最大迭代次数也停止。返回标签向量
5. 溯源找初始传播者：find_source函数，找到初始传播者。遍历图中的节点，如果节点的标签值比所有邻接节点的标签值都大，则认为该节点是初始传播者，将该节点加入到预测感染者列表中并返回。
6. 准确率计算：accuracy函数，这里是调用库一步到位的，f1_score和roc_auc_score函数计算了预测的感染节点与真实感染节点的准确率。
7. 单次模拟：single_simulate函数，上述步骤封装成单次模拟，返回一次模拟的f1和auc的值。在成功单次模拟的基础上，接下来考虑对不同的参数进行优化。
8. 参数优化实验：param_optimize函数，对不同参数（连边概率p、感染概率beta、传播概率alpha、观测时刻t_obs）进行优化实验，控制变量法改变单个参数的值，在每个参数范围内进行maxtimes次模拟（这里取10次模拟的均值作为该参数的f1和auc值,尽量规避实验偶然性），并通过plot函数绘制相关的准确率图表。
连边概率p在范围[0.06, 0.32]内以步长0.02进行计算，每次实验取10次均值。
感染概率beta在范围[0.05, 1.05]内以步长0.05进行计算，每次实验取10次均值。
传播概率alpha在范围[0.1, 1.1]内以步长0.1进行计算，每次实验取10次均值。
观测时刻t_obs在范围[2, 9]内以步长1进行计算，每次实验取10次均值。
9. 主程序：在`__main__`中，设置了实验所需的参数，并依次进行了四次参数优化实验。
10. 由于大部分计算为数值计算，考虑使用并行计算进行优化，这里使用multiprocessing库，获取计算机的CPU数量，最多使用四核心来进行计算，放四个计算任务丢进进程池里，整体速度快了很多。
结果
扫一扫
394
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
