复杂网络建模的实现（哈工大深圳复杂网络建模课程Project）

任务：
1，三张不同的网络：已知某人的名称。已知某人的家乡。已知某人的方言。
分析这三者各自的网络性能（节点度数分布、平均最短路径长度、集聚系数）以及动态行为（在刻意攻击(intentional attack)和随机攻击(random attack)下的鲁棒性）

2，计算三个网络的核数（coreness）

3，画出必要的图片的来显示模拟的结果，以支持观察的结果。

4，（附加内容）设计一个小系统（包括友好的界面以及图形化的展示）来显示整体网络的布局。

节点度数分布
节点i的度数：与该节点所连的边的总数
平均度数：所以节点的度数的平均值
可以直接求

平均最短路径长度：任意两个节点之间的距离的平均值
两节点之间的距离：节点i与节点j之间的边的条数，注意，是最短的边的条数，也就是说两节点之间的距离就是最短的距离（边的条数）
这个需要对图进行分开处理找到的子图

聚集系数：
假设节点i有k(i)条边，也就是说它有k(i)个邻居
定义E(i)：这k(i)个邻居实际有的边的条数
T(i)：这k(i)个邻居可以有的边的条数(k(i)*[K(i)-1]/2)
因此节点i的聚集系数C(i)=E(i)/T(i)

随机攻击：
随机删除节点，查看网络的平均最短路径以及最大连通子图中节点的个数占原先最大连通子图节点个数的比例。

刻意攻击：
每次删除的都是度数最大的节点。评价方式与随机攻击相同。

该实验利用matlab语言进行实现：

• 计算节点的度

function [node_degree]=DegreeDistribution(G)
% 计算每个节点的度

node_degree = sum(G');
figure(1);
bar(node_degree,'r');
title('Node Degree Distribution');
xlabel('Node');
ylabel('Degree');
avg_degree = sum(node_degree)/(length(node_degree));
legend(num2str(avg_degree),'Location','NorthEastOutside');
hold;

figure(2);
md = max(node_degree);
for i=1:md+1
    n_deg(i) = length(find(node_degree == i-1));
end
bar([0:md],n_deg,'r');
title('Node Degree Distribution');
xlabel('Degree');
ylabel('Node');
avg_degree = sum(node_degree)/(length(node_degree));
legend(num2str(avg_degree),'Location','NorthEastOutside');
end

• 计算图的平均最短路径

function [mean_path]=AveragePath( G ,FLAG)
%求一个图的平均最短路径(Average Shortest Path Length)

len=length(G);%节点个数
m = G;
m(m==0)=len*2; %设置为最大距离的两倍，也就是不可达

for i=1:len
    m(i,i)=0;%邻接矩阵的对角线元素均为0
end

if len<=1%少于一个点
    mean_path = 0;
end

for k=1:len%弗洛伊德算法
    for i=1:len
        for j=1:len
            if m(i,j)>m(i,k)+m(k,j)
               m(i,j) = m(i,k)+m(k,j);
            end
        end
    end
end

    
for i=1:len
    m(m==len*2)=0;%不可达点的距离化为0
end
    
   
    each_path = sum(m')/(length(m)-1);%每个点的平均最短路径长度，减一是减去到自己的距离（不存在）
    mean_path = sum(sum(m'))/length(nonzeros(m'));%整个网络的平均最短路径长度
    
    
    if nargin==1%输入的参数为一个，也就是不输入FLAG的值
     bar(each_path,'b');
     title('Average Shortest Path-Length');
     xlabel('Node');
     ylabel('Path length');
     legend(num2str( mean_path),'Location','NorthEastOutside');
    end  

• 找出最大连通子图并计算最短路径

function [max_1,Max_Sub,s,Q,min_path]=MaxSubPath(G,FLAG)
%求最大连通子图的平均最短路径
%C 顶点所分的块数
%con 每行表示一个连通域中的节点标号
%s 各连通域规模（节点数）
%min_path 平均最短路径的长度
%Max_Sub 最大连通子图所在的标号
n=size(G,1);%将矩阵G的行数的值赋给n，也即节点的个数
m=sum(sum(G))/2;
c=1;
q=zeros(n,1);%n*1的全0矩阵
Q=zeros(n,1);
for i=1:n
    if sum(G(i,:))==0  %G整行的和为0，说明i为孤立点
        q(i)=c;
        c=c+1;
    else   
        for j=1:n  
            if G(i,j)==1
                if q(i)==q(j)   %若两者之间有边且标记相同，则属于同一块
                    if q(i)==0
                        q(i)=c;q(j)=c;
                        c=c+1;
                    end
                else    %两者之间有边但标记不同
                    if q(i)==0   %i点尚未获得标记
                        q(i)=q(j);
                    elseif q(j)==0  %j点尚未获得标记
                        q(j)=q(i);
                    else   
                        for k=1:n   %更新标记
                            if q(k)==q(i)
                                q(k)=q(j);
                            end
                        end
                    end
                end
            end
        end
    end
end

%重新编号，排除空的块号
C=0;
for k=1:c-1
    if ~isempty(find(q==k, 1))%利用find函数找到q==k的下标
        C=C+1;
        Q(q==k)=C;%存储各个节点所在的块号
        s(C,1)=length(find(q==k)); 
    end
end

sort(s,1,'descend');

Max_Sub=find(s==max(s));%最大连通子图所在的标号
len_Q = length(Q);


M = zeros(len_Q,len_Q)%创建一个新的矩阵存储最大连通子图的邻接矩阵
for i=1:len_Q
  if Q(i)==Max_Sub
      for j=1:len_Q
          M(i,j)=G(i,j);
      end
  end
end

max_1 = s(1);%找到最大连通子图的节点数
   for i=1:length(s)
       if s(i)>max_1
           max_1 = s(i);
       end
   end
   
A1=sum(abs(M'));%删除所有为0的行
index=find(A1==0);
M(index,:)=[];

A2=sum(abs(M));%删除所有为0的列
index=find(A2==0);
M(:,index)=[];

if nargin==1
    AveragePath(M);%计算最大连通子图的平均最短路径
else
    min_path = AveragePath(M,FLAG);%计算最大连通子图的平均最短路径
end  

• 计算各节点的聚集系数

function [C]=Coefficient(G)
%计算每个节点的聚集系数(Cluster Coefficient)

K = sum(G');%每个节点的度
C = zeros(1,length(G));%记录每个节点的聚集系数

for i=1:length(G)
    a = find(G(i,:)>0);
    E=0;
    for j=1:length(a)
      for k=j+1:length(a)
          if G(a(j),a(k))==1
            E = E+1;%邻居已有的边的条数
          end
      end
    end
    if K(i)==1||K(i)==0
        C(i)=0;
    else
        C(i) = 2*E/(K(i)*(K(i)-1));
    end
end


mean_c = sum(C)/length(C);
bar(C,'g');
title('Clustering Coefficient');
xlabel('Node');
ylabel('Coefficient');
legend(num2str(mean_c),'Location','NorthEastOutside');
end

• 计算网络的核数

function [CORENESS]=Coreness(G)
%计算每个节点以及整个网络的核数(Coreness)
CORENESS = zeros(1,length(G));%1*n的全0矩阵，即所有节点的核数均为0
nodes_degree = sum(G');%每个点的度
max = length(G)+ 1;%max即节点所能达到的最大度数+1，也即不能达到的度数
pre_min = 0;%负责记录上一轮中的最小的度数
counted_num=0;%负责记录查找并计算过的节点的数

for i=1:length(nodes_degree)
   m = min(nodes_degree);%返回度数的最小值，是一个数值
   min_node = find(nodes_degree==m);%是一个列向量,返回最小度数的节点所在的列的下标，也即度数最小的节点序号
   counted_num = counted_num + length(min_node);
   if pre_min>m
   CORENESS(min_node) = pre_min;
   else
       CORENESS(min_node) = m;
   end
   
   if length(min_node)>1     %一个以上的最小度数的点
      for j=1:length(min_node)
          x = min_node(j);
          for k=1:length(G)
              if G(k,x)==1%节点k与最小度数节点相连            
                  if CORENESS(k)==0%且该节点尚未在计算核数的范围内
                    nodes_degree(k) = nodes_degree(k)-1;%更新度数
                  end
              end
          end
      end
   else     %只有一个最小度数的点
       for j=1:length(G)
          if G(j,min_node)==1;
              if CORENESS(j)==0
                 nodes_degree(j) = nodes_degree(j)-1;%更新度数
              end
          end
       end
   end
   
   nodes_degree(min_node) = max;%将这些原先度数最小的点的度化为max，目的是为了下一轮找最小度数节点不会再找到它们
   pre_min = m;%记录上一轮中的最小度数
   if counted_num>=length(nodes_degree)%查找完毕，提前跳出整个for循环
       break;
   end
end


bar(CORENESS,'m');
title('Node Coreness');
xlabel('Node');
ylabel('Coreness');


max_coreness = CORENESS(1);%找到最大的核数作为整个网络的核数
for i=1:length(CORENESS)
    if CORENESS(i)>max_coreness
        max_coreness = CORENESS(i);
    end
end

legend(num2str(max_coreness),'Location','NorthEastOutside');

end

• 随机攻击

function [mean_path]=RandomAttack1( G )
%对网络进行随机攻击，查看其对平均最短路径的影响(RandomAttack)
mean_path = zeros(1,length(G));
[max,~,~,~,~] = MaxSubPath(G,2);
MaxSubNum = max;
max_sub_num = zeros(1,length(G));
temp = G;


for i=1:length(G)-1
   temp(1,:)=[];%去掉第1行
   temp(:,1)=[];%去掉第1列
   [~,~,~,~,min_path] = MaxSubPath(temp,2);
   mean_path(i) = min_path;
   [max_1,~,~,~,~] = MaxSubPath(temp,2);
   max_sub_num(i) = max_1/MaxSubNum;
end

figure(1);
x = (1:length(G));
plot(x,mean_path,'bx');
title('Robustness Against Random Attack');
xlabel('Node');
legend('Average Shortest Path of Largest SubGraph','Location','NorthEastOutside');
hold on;
figure(2);
plot(x,max_sub_num,'ro')
title('Robustness Against Random Attack');
xlabel('Node');
legend('Percentage of Largest SubGraph','Location','NorthEastOutside');
end

• 刻意攻击

function [keep_path]=IntentionalAttack1( G)
%对网络进行随机攻击，查看其对平均最短路径的影响(IntentionalAttack)
degree = sum(G');
mean_path = zeros(1,length(G));
[max_2,~,~,~,~] = MaxSubPath(G,2);
MaxSubNum = max_2;
max_sub_num = zeros(1,length(G));
temp = G;

for i=1:length(G)-1
  degree = sum(temp');
  j = find(degree==max(degree));
  if length(j)>1%不止一个的话就取第一个
      j = j(1);
  end
  temp(j,:)=[];
  temp(:,j)=[]; 
   
   [~,~,~,~,min_path] = MaxSubPath(temp,2);
   mean_path(i) = min_path;
   [max_1,~,~,~,~] = MaxSubPath(temp,2);
   max_sub_num(i) = max_1/MaxSubNum;
   
end


figure(1);
x = (1:length(G));
plot(x,mean_path,'bx');
title('Robustness Against Intentional Attack');
xlabel('Node');
legend('Average Shortest Path of Largest SubGraph','Location','NorthEastOutside');
hold on;
figure(2);
plot(x,max_sub_num,'ro')
title('Robustness Against Intentional Attack');
xlabel('Node');
legend('Percentage of Largest SubGraph','Location','NorthEastOutside');
end