scale_free.m

最短距离的求解

% 2006-12-20
% 无标度网络仿真
tic;

% 四个输入变量
I = 30 ;    % 表示现实的次数,I要大于或者等于3,才能对节点的度数用邻接矩阵进行统计
N = 1000 ;  % 表示网络的节点的个数
m0 = 3 ;    % 表示网络的初始节点个数
m = 3 ;     % 表示新点与旧网络连边的数目

% 输出变量realization_of_distribution
realization_of_distribution = sparse( I , N ) ;  % 矩阵的每行 存储 度分布的一个现实

for J = 1 : I   % 对I次现实做平均，然后用这个平均值近似网络的度分布
format long;

adjacent_matrix = sparse( m0 , m0 ) ;  % 用sparse表示邻接矩阵adjacent_matrix，极大地释放内存
for i = 1 : m0
    for j = 1 : m0
        if j ~= i
            adjacent_matrix( i , j ) = 1 ;
        end
    end
end
adjacent_matrix = sparse( adjacent_matrix ) ;  % 初始网络，这里利用sparse把内存得以释放

% 修改cumsum,目的是使得网络的顶点的个数能够超过10万,而不超过内存空间的限制
node_degree = sparse( 1 , m0 ) ;  % node_degree表示各个节点的度数
for p = 1 : m0
    node_degree( p ) = sum( adjacent_matrix( 1 : m0 , p ) ) ;
end

% 每次加入一个新点，新点和老点之间按照择优概率进行连接，值得注意的是，每次新点加入之后，
% 网络的择优概率都发生变化，每一次循环都是一个相对独立的整体，要把流程进行分割处理
for iteration = m0 + 1 : N
    [ J , iteration ]  % 控制现实和迭代的次数    

    % node_degree是每次循环所需要的唯一输出变量    
    % 为每次迭代配置相关的变量
    total_degree = 2 * m * ( iteration - 4 ) + 6 ;  % 迭代之前的网络各个节点的度数之和
    degree_frequency = node_degree / total_degree ;  % 每个节点的度数的频数，这是新点连边的择优概率
    cum_distribution = cumsum( degree_frequency ) ;  % cum_distribution把区间(0,1]分成若干个小区间，从而对这些个小区间进行投点实验
    
    % 下面把 r1 和 (0,1] 内的各个小区间的端点进行比较，落在第i小区间，就意味着和第i个节点相连边
    choose = zeros( 1 , m ) ;  % choose存放的是和新点相连接的三个老点
    r1 = rand(1) ;    % rand 的范围为(0,1)
    choose(1) = min( find( ( cum_distribution >= r1) == 1 ) ) ;

    % 选出第二个和新点相连接的顶点，注意这两个点是不同的，目的是避免重复边的出现
    r2 = rand(1) ;
    choose(2) = min( find( ( cum_distribution >= r2) == 1 ) ) ;
    while choose(2) == choose(1)
        r2 = rand(1) ;
        choose(2) = min( find( ( cum_distribution >= r2) == 1 ) ) ;
    end

    % 选出第三个和新点相连接的顶点，注意这三个点是不同的，目的是避免重复边的出现
    r3 = rand(1) ;
    choose(3) = min( find( ( cum_distribution >= r3) == 1 ) ) ;
    while ( choose(3) == choose(1) ) | ( choose(3) == choose(2) )
        r3 = rand(1) ;
        choose(3) = min( find( ( cum_distribution >= r3) == 1 ) ) ;
    end

    % 把新点加入网络后，对邻接矩阵进行相应的改变！
    % 这是在一次循环下生成的新的邻接矩阵，下一次循环就是在这个邻接矩阵的基础上进行的！
    for k = 1 : m
        adjacent_matrix( iteration , choose(k) ) = 1 ;
        adjacent_matrix( choose(k) , iteration ) = 1 ;
    end

    for p = 1 : iteration
        node_degree(p) = sum( adjacent_matrix( 1 :iteration , p ) ) ;  % 这个循环的目的是重新给各个节点的度赋值
    end    
end  % //for iteration = 4 : N  一次最外层循环的结束

number_of_nodes_with_equal_degree = zeros( 1 , N ) ;    % 存储度相同的顶点的个数

for i = 1 : N
    difference = node_degree - i * ones( 1 , N ) ;
    number_of_nodes_with_equal_degree(i) = length( find( difference == 0 ) ) ;  % 度为i的节点的个数
end
realization_of_distribution( J , : ) = number_of_nodes_with_equal_degree ;
end  % //for J = 1 : I  外层循环的中止

aaa = cumsum( realization_of_distribution ) ;
bb1 = aaa( I , : ) ;   % 譬如，度为3的节点的个数，由于度数为1,2的节点的个数为0,故可以从度数为3的节点个数开始计算
bb2 = bb1( m : N ) ;
bbb = bb2 / ( I * N ) ;  % 譬如，度为3的节点的个数在网络中的比例
K = m : N ;
loglog( K , bbb , '*' ) ; % 注意，作图的时候，一定要做散点图
axis([1 N 0.0000001 0.9]);
hold on;
y1 = 2 * m^2 * K .^ ( -3 ) ;
loglog( K , y1 , 'r' ) ;  % 与平均场结果进行比较 p(k)=2*m^2*k^(-3)
  
% 第四步：：全部工作结束
toc;  % 计算程序运行需要的时间