论文已中,把刚开始写的代码放一下(最终版就不放了),版权所有,禁止抄袭,转载请注明出处。
目标函数
网络效率为 η = 1 N ( N − 1 ) ∑ v i ≠ v j ∈ V 1 d i j \eta =\frac{1}{{N(N - 1)}}\sum\limits_{v_i \ne v_j \in V}\frac{1}{{d_{ij}}} η=N(N−1)1vi=vj∈V∑dij1
则衡量网络毁伤的目标函数即网络下降效率为
R
(
s
)
=
(
1
−
η
(
G
′
)
/
η
(
G
)
)
×
100
%
R(s) = (1-\eta (G')/\eta(G))\times100\%
R(s)=(1−η(G′)/η(G))×100%
下面对比了基于拓扑势的CELF算法(TPCELF),基于重要节点的贪心算法(GABIN)和基于度值(Degree)、接近度(Closeness)和介数(Betweenness)排序的算法,Matlab程序是初步尝试,犹豫在计算介数的时候总存在一些时间消耗过大的问题,按照常用的计算介数方法编程消耗太大,对大型网络等候时间过长,Lev的工具箱暂时不好用,总是报错。这里改用了matlab_bgl工具箱,具体的地址见这里,使用的时候需要把private文件夹重命名一下在添加到路径。Matlab的实验比较简单,基于拓扑势的方法也没有加CELF。
后来发现python的好几个包在这块计算上做的比较好,改用基于python编程。
代码总体是完整的,只是没有放生成仿真数据的的内容,在另一个文档里,想比读者肯定能够自己实现的,非常easy,网上也有很多现成的数据集,改一下路径就能跑~
Matlab代码
close all
clear
tic
%% 生成有无标度网络
%方式一 自己生成
% BA=Powerlaw_with_Expected_Exponent(300,6,2,1);%这个函数生成的幂律分布的网络会有孤立点存在
%方式二 用snap网络生成的数据集,如果是批量的,加一个循环就可以
m = 1;
filename = strcat('D:\test_300_1.5_',num2str(m),'.txt');
d=textread(filename,'','headerlines',3);%跳过第3行开始读数据
s=[];t=[];
for i=1:length(d)
s=[s,d(i,1)+1];
t=[t,d(i,2)+1];
end
g = graph(s,t);
BA = adjacency(g); %下文介数是要用稀疏矩阵的
A = BA;
A_1 = BA; %用于比较,保存初始矩阵,liu的方法
G =graph(BA);
% BA = adjacency(g); %得到邻接矩阵
D = distances(G); %得到距离矩阵
plot(G) %展示图片是怎样的
%% 初始化相应的值
N=height(G.Nodes); %节点数
DD = D+1e6*eye(N); %在对角线上为了排除自己和自己的距离,令其为一个很大的数
Initial_efficiency = sum(sum(1./DD)); %初始的网络效率
S=[]; S_1=[]; %两个种子集
sigma = 0.9428; %影响力因子,此处设置为两跳的最小值
alpha = 8; %Liu论文的节点初始集合放大系数
%% 想要打击K个节点
K=5;
for k=K: -1 :1
%计算节点的拓扑势
for i=1:N
phi(i)=sum(exp(-(D(i,:)/(sigma)).^2));
end
[BB,I] = sort( phi,'descend'); %对拓扑势进行排序
[B,E] = betweenness_centrality(A_1); %计算介数
Du = Degree_undirect(A_1); %计算度数
De = Du+max(B)/max(Du).*B; %计算节点重要度
[BB_1,I_1] = sort( De ,'descend'); %对重要度进行排序
%这里取用最小alpha=10,相对固定住
for i = 1:max(k,10)
a=A;
a(I(i),:)=0;
a(:,I(i))=0;
g = graph(a);
D = distances(g);
DD = D+1e6*eye(N);
descend_efficiency(i) = 1-sum(sum(1./DD))/Initial_efficiency;
end
for i=1:alpha*k
a_1=A_1;
a_1(I(i),:)=0;
a_1(:,I(i))=0;
g_1 = graph(a_1);
D_1 = distances(g_1);
DD_1 = D_1+1e6*eye(N);
descend_efficiency_1(i) = 1-sum(sum(1./DD_1))/Initial_efficiency;
end
ddd(K+1-k)=max(descend_efficiency);
index=I(find(descend_efficiency == max(descend_efficiency),1));%找到最大衰减性能对应的节点的位置
S = [S,index]; %把节点加入种子集合
A(index,:)=0; %删除被打击的节点
A(:,index)=0;
ddd_1(K+1-k)=max(descend_efficiency_1);
index_1=I_1(find(descend_efficiency_1 == max(descend_efficiency_1),1));%找到最大衰减性能对应的节点的位置
S_1 = [S_1,index_1]; %把节点加入种子集合
A_1(index_1,:)=0; %删除被打击的节点
A_1(:,index_1)=0;
end
toc
figure
plot(ddd,'-*');hold on ;
plot(ddd_1)
legend('我','Liu')
Python代码
import networkx as nx
import matplotlib.pyplot as plt
import time
import numpy as np
import copy
import pandas as pd
import os
'''
=====计算增益的等价说明======
在计算每个节点删除后的增益时,采用了选择最大下降效率的方法,其实这和最大增益是等价的
因为最大下降效率 - 上一状态网络的下降效率 = 最大增益,被减数在每轮中是一个定值
'''
def txt_read(data):
'''
读取txt文档,生成网络
:param data:
:return:
'''
return nx.read_edgelist(data, create_using=nx.Graph())
def effectiveness(path):
'''
读取最短路径,计算网络的效率count
:param path:
:return:
'''
count = 0
for x in path.keys():
for y in path[x].keys():
try:
count += 1/path[x][y]
except:
pass
return count
def CELF_algrithm(path, Use):
'''
CELF算法封装
传入总体参数 g N
Initial_efficiency
:param path: 初始的最短路径
:param Use: bool 是否用celf思想
:return: S_1 decrease_efficiency_1 time1
'''
# ##### #####开始进行网络毁伤实验 Part1 ##### #####
start1 = time.time()
decrease_efficiency_1 = []
S_1 = []
k = 0
while k < K:
phi = []
for i in range(N): # 这样的用法是为了让节点按照顺序排列
temp = 0
for j in g.nodes():
try:
temp += np.exp(-(path[j][str(i)] / (sigma)) ** 2)
except KeyError:
# 有的节点之间没有最短路径,会报错,采用异常处理模式,这样就不影响结果
pass
phi.append(temp)
idx = sorted(enumerate(phi), key=lambda x: x[1], reverse=True) # 降序排列,得到id的元组
efficiency = []
for p in range(max(K - k, flex_alpha)):
a = copy.deepcopy(G_1)
a.remove_node(str(idx[p][0])) # 删除节点以及对应的边
a.add_node(str(idx[p][0])) # 添加被删除的节点
# 得到网络的最短路径,保存为两级的字典格式
path = dict(nx.shortest_path_length(a))
# 计算网络下降效率
efficiency.append(1 - effectiveness(path) / Initial_efficiency)
id_top3 = pd.Series(efficiency).sort_values(ascending=False).index[:3] # 获取效率排名前三的效率值的id
decrease_efficiency_1.append(efficiency[id_top3[0]]) # 选择效率下降最大的节点,也可用max(efficiency)
index = str(idx[efficiency.index(efficiency[id_top3[0]])][0]) # 下标是字符串类型的,效率下降对应的p值,再对应返回过去idx
S_1.append(index)
G_1.remove_node(index) # 删除节点以及对应的边
G_1.add_node(index) # 添加被删除的节点
# plt.figure()
# nx.draw_networkx(G_1)
if Use == True:
# 采用 CELF 算法思想:B的新一轮的增益要大于C在本轮的增益
# 判断一下是否已经到达零界位置,否则在k刚好满足要求的时候会多运行一步
if k + 1 == K:
break
# 计算C在本轮的增益
if k == 0: # 如果是打击第一个节点
delta_C = efficiency[id_top3[2]]
# print('delta_C=', delta_C)
else:
delta_C = efficiency[id_top3[2]] - decrease_efficiency_1[k - 1]
# print('delta_C=',delta_C)
index_B_new = id_top3[1]
a = copy.deepcopy(G_1)
a.remove_node(str(idx[index_B_new][0])) # 删除节点以及对应的边
a.add_node(str(idx[index_B_new][0])) # 添加被删除的节点
# 得到网络的最短路径,保存为两级的字典格式
path = dict(nx.shortest_path_length(a))
# 计算网络下降效率增益
delta_B_new = (1 - effectiveness(path) / Initial_efficiency) - decrease_efficiency_1[k]
if delta_C < delta_B_new:
# 如果新的一轮B的增益大于本轮C的增益,直接用B,然后删除节点
decrease_efficiency_1.append(1 - effectiveness(path) / Initial_efficiency)
S_1.append(str(idx[index_B_new][0]))
G_1.remove_node(str(idx[index_B_new][0])) # 删除节点以及对应的边
G_1.add_node(str(idx[index_B_new][0])) # 添加被删除的节点
k += 2
# plt.figure()
# nx.draw_networkx(G_1)
else: #如果B的新增益比原来的C要小
k += 1
else: #如果不同celf思想
k += 1
end1 = time.time()
time1 = end1 - start1
return S_1, decrease_efficiency_1, time1
def GABIN_algrithm():
'''
:return:
'''
# ##### #####开始进行网络毁伤实验 Part2 ##### #####
start2 = time.time()
decrease_efficiency_2 = []
S_2 = []
k = 0
while k < K:
# 得到网络的介数,存储是字典
B = nx.betweenness_centrality(G_2, normalized=False)
# 得到的度值,存储是字典
Du = nx.degree(G_2)
# 计算最大介数和最大度之间的比值,重要度
De = max(B.values()) / max(Du, key=lambda x: x[-1])[-1]
# 计算节点重要度
Importance_D_B = {}
for key in B.keys():
Importance_D_B[key] = Du[key] + De * B[key]
# [Importance_D_B[key]] = [Du[key] + De * B.get(key) for key in B.keys()] # 有问题,这样就把原来的序号给搞没了
idx = sorted(Importance_D_B.items(), key=lambda x: x[1], reverse=True) # 降序排列
efficiency = []
for p in range(alpha * (K - k)):
a = copy.deepcopy(G_2)
a.remove_node(str(idx[p][0])) # 删除节点以及对应的边
a.add_node(str(idx[p][0])) # 添加被删除的节点
# 得到网络的最短路径,保存为两级的字典格式
path = dict(nx.shortest_path_length(a))
# 计算网络下降效率
efficiency.append(1 - effectiveness(path) / Initial_efficiency)
decrease_efficiency_2.append(max(efficiency)) # 选择效率下降最大的节点
index = str(idx[efficiency.index(max(efficiency))][0]) # 下标是字符串类型的,效率下降对应的p值,再对应返回过去idx
S_2.append(index)
G_2.remove_node(index) # 删除节点以及对应的边
G_2.add_node(index) # 添加被删除的节点
k += 1
end2 = time.time()
time2 = end2 - start2
return S_2, decrease_efficiency_2, time2
def Greedy_algrithm():
'''
贪婪算法
:param g:
:return:
'''
# ##### #####开始进行网络毁伤实验 Part3 ##### #####
start3 = time.time()
decrease_efficiency_3 = []
S_3 = []
k = 0
while k < K:
efficiency = []
for p in G_3.nodes():
a = copy.deepcopy(G_3)
a.remove_node(p) # 删除节点以及对应的边
a.add_node(p) # 添加被删除的节点
# 得到网络的最短路径,保存为两级的字典格式
path = dict(nx.shortest_path_length(a))
# 计算网络下降效率,保存成节点编号和效率值的tuple
efficiency.append((p, 1 - effectiveness(path) / Initial_efficiency))
decrease_efficiency_3.append(max(efficiency,key= lambda x:x[-1])[-1]) # 选择效率下降最大的节点
index = max(efficiency,key= lambda x:x[-1])[0] # 下标是字符串类型的
S_3.append(index)
G_3.remove_node(index) # 删除节点以及对应的边
G_3.add_node(index) # 添加被删除的节点
k += 1
end3 = time.time()
time3 = end3 - start3
return S_3, decrease_efficiency_3, time3
def Du_algrithm():
'''
:return:
'''
start4 = time.time()
G_4 = copy.deepcopy(g)
Du = nx.degree(G_4)
temp = sorted(Du, key=lambda x: x[1], reverse=True)[0:K] # 元组降序排列
S_4 = [temp[key][0] for key in range(K)] # 得到种子集
decrease_efficiency_4 = []
for p in S_4:
G_4.remove_node(p)
G_4.add_node(p)
path = dict(nx.shortest_path_length(G_4))
decrease_efficiency_4.append(1 - effectiveness(path) / Initial_efficiency)
end4 = time.time()
def Closness_algrithm():
'''
:return:
'''
start5 = time.time()
G_5 = copy.deepcopy(g)
Closeness = nx.closeness_centrality(G_5)
temp = sorted(Closeness.items(), key=lambda x: x[1], reverse=True)[0:K] # 字典降序排列
S_5 = [temp[key][0] for key in range(K)] # 得到种子集
decrease_efficiency_5 = []
for p in S_5:
G_5.remove_node(p)
G_5.add_node(p)
path = dict(nx.shortest_path_length(G_5))
decrease_efficiency_5.append(1 - effectiveness(path) / Initial_efficiency)
end5 = time.time()
return S_5, decrease_efficiency_5
def Betweenness_algrithm():
'''
:return:
'''
start6 =time.time()
G_6 = copy.deepcopy(g)
Betweenness = nx.betweenness_centrality(G_6)
temp = sorted(Betweenness.items(), key=lambda x: x[1], reverse=True)[0:K] # 字典降序排列
S_6 = [temp[key][0] for key in range(K)] # 得到种子集
decrease_efficiency_6 = []
for p in S_6:
G_6.remove_node(p)
G_6.add_node(p)
path = dict(nx.shortest_path_length(G_6))
decrease_efficiency_6.append(1 - effectiveness(path) / Initial_efficiency)
end6 = time.time()
return S_6, decrease_efficiency_6
if __name__ == '__main__':
'''
实验1获取数据 设定gamma=2 K=5 alpha=5 flex_alpha=10
####################################################
'''
Num_of_Nodes = [50, 60, 70, 80, 90, 100, 200, 300, 400, 500, 600, 700]
Gamma = 2
fsave = open('result1.txt', mode='w+')
for num in Num_of_Nodes:
filename = 'D:\\PY\\TEST1\\test_' + str(num) + '_2.txt'
g = txt_read(filename)
'''
####################################################
'''
# # 得到网络的最短路径,保存为两级的字典格式
path = dict(nx.shortest_path_length(g))
# # 设置备份网络并展示原始网络,不能简单传递变量赋值,这是链接到内存地址的
G_1 = copy.deepcopy(g)
G_2 = copy.deepcopy(g)
G_3 = copy.deepcopy(g)
'''
##### ##### ##### ##### ##### ##### ##### #####
初始化参数设置
'''
sigma = 0.9428
# 节点个数
N = g.number_of_nodes()
# 需要毁伤的节点个数
K = 5
# 种子候选点的放大系数
alpha = 5
# 拓扑势方法单一选择的最大的值
flex_alpha = 10
# 初始的网络效率
Initial_efficiency = effectiveness(path)
'''
##### ##### ##### ##### ##### ##### ##### #####
'''
# 输出部分的内容
out1 = CELF_algrithm(path, True)
out2 = GABIN_algrithm()
out3 = Greedy_algrithm()
fsave.write(str(num)+', '+str(out1[2])+', '+str(out2[2])+', '+str(out3[2])+'\n')
fsave.close()
start = time.time()
'''
实验2 获取数据 设定gamma=2.5 N=300 K=1-10 alpha=5 flex_alpha=10
####################################################
'''
# 设置K
K = 10
fsave = open('K_=_10.txt', 'w+')
for i in range(1, 101):
filename = 'D:\\PY\\TEST2\\test_300_' + str(2.5) + '_' + str(i) + '.txt'
g = txt_read(filename)
# 得到网络的最短路径,保存为两级的字典格式
path = dict(nx.shortest_path_length(g))
# 设置备份网络并展示原始网络,不能简单传递变量赋值,这是链接到内存地址的
G_1 = copy.deepcopy(g)
G_2 = copy.deepcopy(g)
'''
##### ##### ##### ##### ##### ##### ##### #####
初始化参数设置
'''
sigma = 0.9428
# 节点个数
N = g.number_of_nodes()
# 种子候选点的放大系数
alpha = 5
# 拓扑势方法单一选择的最大的值
flex_alpha = 10
# 初始的网络效率
Initial_efficiency = effectiveness(path)
'''
##### ##### ##### ##### ##### ##### ##### #####
'''
out1 = CELF_algrithm(path, True)
out2 = GABIN_algrithm()
out3 = Du_algrithm()
out4 = Closness_algrithm()
out5 = Betweenness_algrithm()
fsave.write(str(out1) + '\n')
fsave.write(str(out2) + '\n')
fsave.write(str(out3) + '\n')
fsave.write(str(out4) + '\n')
fsave.write(str(out5) + '\n')
fsave.write('\n')
fsave.close()
end = time.time()
print(end-start)