% 2006-12-20
% 無標度網(wǎng)絡(luò)仿真
tic;

% 四個輸入變量
I = 30 ;    % 表示現(xiàn)實的次數(shù),I要大于或者等于3,才能對節(jié)點的度數(shù)用鄰接矩陣進行統(tǒng)計
N = 1000 ;  % 表示網(wǎng)絡(luò)的節(jié)點的個數(shù)
m0 = 3 ;    % 表示網(wǎng)絡(luò)的初始節(jié)點個數(shù)
m = 3 ;     % 表示新點與舊網(wǎng)絡(luò)連邊的數(shù)目

% 輸出變量realization_of_distribution
realization_of_distribution = sparse( I , N ) ;  % 矩陣的每行 存儲 度分布的一個現(xiàn)實

for J = 1 : I   % 對I次現(xiàn)實做平均，然后用這個平均值近似網(wǎng)絡(luò)的度分布
format long;

adjacent_matrix = sparse( m0 , m0 ) ;  % 用sparse表示鄰接矩陣adjacent_matrix，極大地釋放內(nèi)存
for i = 1 : m0
    for j = 1 : m0
        if j ~= i
            adjacent_matrix( i , j ) = 1 ;
        end
    end
end
adjacent_matrix = sparse( adjacent_matrix ) ;  % 初始網(wǎng)絡(luò)，這里利用sparse把內(nèi)存得以釋放

% 修改cumsum,目的是使得網(wǎng)絡(luò)的頂點的個數(shù)能夠超過10萬,而不超過內(nèi)存空間的限制
node_degree = sparse( 1 , m0 ) ;  % node_degree表示各個節(jié)點的度數(shù)
for p = 1 : m0
    node_degree( p ) = sum( adjacent_matrix( 1 : m0 , p ) ) ;
end

% 每次加入一個新點，新點和老點之間按照擇優(yōu)概率進行連接，值得注意的是，每次新點加入之后，
% 網(wǎng)絡(luò)的擇優(yōu)概率都發(fā)生變化，每一次循環(huán)都是一個相對獨立的整體，要把流程進行分割處理
for iteration = m0 + 1 : N
    [ J , iteration ]  % 控制現(xiàn)實和迭代的次數(shù)    

    % node_degree是每次循環(huán)所需要的唯一輸出變量    
    % 為每次迭代配置相關(guān)的變量
    total_degree = 2 * m * ( iteration - 4 ) + 6 ;  % 迭代之前的網(wǎng)絡(luò)各個節(jié)點的度數(shù)之和
    degree_frequency = node_degree / total_degree ;  % 每個節(jié)點的度數(shù)的頻數(shù)，這是新點連邊的擇優(yōu)概率
    cum_distribution = cumsum( degree_frequency ) ;  % cum_distribution把區(qū)間(0,1]分成若干個小區(qū)間，從而對這些個小區(qū)間進行投點實驗
    
    % 下面把 r1 和 (0,1] 內(nèi)的各個小區(qū)間的端點進行比較，落在第i小區(qū)間，就意味著和第i個節(jié)點相連邊
    choose = zeros( 1 , m ) ;  % choose存放的是和新點相連接的三個老點
    r1 = rand(1) ;    % rand 的范圍為(0,1)
    choose(1) = min( find( ( cum_distribution >= r1) == 1 ) ) ;

    % 選出第二個和新點相連接的頂點，注意這兩個點是不同的，目的是避免重復(fù)邊的出現(xiàn)
    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

    % 選出第三個和新點相連接的頂點，注意這三個點是不同的，目的是避免重復(fù)邊的出現(xiàn)
    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

    % 把新點加入網(wǎng)絡(luò)后，對鄰接矩陣進行相應(yīng)的改變！
    % 這是在一次循環(huán)下生成的新的鄰接矩陣，下一次循環(huán)就是在這個鄰接矩陣的基礎(chǔ)上進行的！
    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 ) ) ;  % 這個循環(huán)的目的是重新給各個節(jié)點的度賦值
    end    
end  % //for iteration = 4 : N  一次最外層循環(huán)的結(jié)束

number_of_nodes_with_equal_degree = zeros( 1 , N ) ;    % 存儲度相同的頂點的個數(shù)

for i = 1 : N
    difference = node_degree - i * ones( 1 , N ) ;
    number_of_nodes_with_equal_degree(i) = length( find( difference == 0 ) ) ;  % 度為i的節(jié)點的個數(shù)
end
realization_of_distribution( J , : ) = number_of_nodes_with_equal_degree ;
end  % //for J = 1 : I  外層循環(huán)的中止

aaa = cumsum( realization_of_distribution ) ;
bb1 = aaa( I , : ) ;   % 譬如，度為3的節(jié)點的個數(shù)，由于度數(shù)為1,2的節(jié)點的個數(shù)為0,故可以從度數(shù)為3的節(jié)點個數(shù)開始計算
bb2 = bb1( m : N ) ;
bbb = bb2 / ( I * N ) ;  % 譬如，度為3的節(jié)點的個數(shù)在網(wǎng)絡(luò)中的比例
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' ) ;  % 與平均場結(jié)果進行比較 p(k)=2*m^2*k^(-3)
  
% 第四步：：全部工作結(jié)束
toc;  % 計算程序運行需要的時間