%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 生成scale-free網(wǎng)絡(luò)，目的是在這個(gè)網(wǎng)絡(luò)結(jié)構(gòu)上展開相關(guān)的討論！
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tic;
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% 開發(fā)人：： 于灝
% 單位：： 東北大學(xué)信息科學(xué)與工程學(xué)院
% 開發(fā)日期：： 2006年9月11日
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%**************************************************************************
%%% 網(wǎng)絡(luò)的全局變量
%**************************************************************************
%%%%%% 四個(gè)輸入變量
I=10; %% 表示現(xiàn)實(shí)的次數(shù),I要大于或者等于3,才能對(duì)節(jié)點(diǎn)的度數(shù)用鄰接矩陣進(jìn)行統(tǒng)計(jì)
N=20; %% 表示網(wǎng)絡(luò)的節(jié)點(diǎn)的個(gè)數(shù)
m0=3; %% 表示網(wǎng)絡(luò)的初始節(jié)點(diǎn)個(gè)數(shù)
m=3; %% 表示新點(diǎn)與舊網(wǎng)絡(luò)連邊的數(shù)目
%%%%%% 只有一個(gè)輸出變量realization_of_distribution
realization_of_distribution=sparse(I,N); % 矩陣的每行存儲(chǔ)度分布的一個(gè)現(xiàn)實(shí)
%**************************************************************************
%**************************************************************************
%**************************************************************************
%**************************************************************************
for J=1:I % 對(duì)I次現(xiàn)實(shí)做平均，然后用這個(gè)平均值近似網(wǎng)絡(luò)的度分布
format long;
adjacent_matrix=sparse(m0,m0); % 用sparse表示鄰接矩陣adjacent_matrix，極大地釋放內(nèi)存
for i=1:m0
for j=1:3
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ò)的頂點(diǎn)的個(gè)數(shù)能夠超過(guò)10萬(wàn),而不超過(guò)內(nèi)存空間的限制
node_degree=sparse(1,m0); % node_degree表示各個(gè)節(jié)點(diǎn)的度數(shù)
for p=1:m0
%last_element=sparse(m0,1);
%last_element=cumsum(adjacent_matrix(1:m0,p));
%node_degree(p)=last_element(m0);
node_degree(p)=sum(adjacent_matrix(1:m0,p));
end
%%%%%%%%%%%%%%%%%%%%%%%%%
%**************************************************************************
%**************************************************************************
%**************************************************************************
% 每次加入一個(gè)新點(diǎn)，新點(diǎn)和老點(diǎn)之間按照擇優(yōu)概率進(jìn)行連接，
% 值得注意的是，每次新點(diǎn)加入之后，網(wǎng)絡(luò)的擇優(yōu)概率都發(fā)生變化，
% 每一次循環(huán)都是一個(gè)相對(duì)獨(dú)立的整體，要把流程進(jìn)行分割處理
for iteration=4:N
[J,iteration] % 控制現(xiàn)實(shí)和迭代的次數(shù) 
%**************************************************************************
%%%% node_degree是每次循環(huán)所需要的唯一輸出變量 
%**************************************************************************
%%% 為每次迭代配置相關(guān)的變量
total_degree=2*m*(iteration-4)+6; %%% 迭代之前的網(wǎng)絡(luò)各個(gè)節(jié)點(diǎn)的度數(shù)之和
degree_frequency=node_degree/total_degree; %%% 每個(gè)節(jié)點(diǎn)的度數(shù)的頻數(shù)，這是新點(diǎn)連邊的擇優(yōu)概率
cum_distribution=cumsum(degree_frequency); %%% cum_distribution把區(qū)間 [0,1] 分成若干個(gè)小區(qū)間，從而對(duì)這些個(gè)小區(qū)間進(jìn)行投點(diǎn)實(shí)驗(yàn)
interval=cum_distribution(1:(iteration-1)); %%% 這是小區(qū)間的端點(diǎn)，是cum_distribution的前 iteration－1 個(gè)元素
%**************************************************************************
%%% 下面把 r1 和 [0,1] 內(nèi)的各個(gè)小區(qū)間的端點(diǎn)進(jìn)行比較，落在第 i 小區(qū)間，就意味著和第 i 個(gè)節(jié)點(diǎn)相連邊 %%% 
choose=zeros(1,m); %%% choose存放的是和新點(diǎn)相連接的三個(gè)老點(diǎn) 
%r2=rand(1); %||| 任意選擇的一個(gè) [0,1] 隨機(jī)數(shù)
%r3=rand(1); %|||
%**************************************************************************
%**************************************************************************
%%% 選出第一個(gè)和新點(diǎn)相連接的頂點(diǎn)
r1=rand(1); 
if r1 choose(1)=1;
elseif r1>=interval(iteration-2)
choose(1)=iteration-1;
elseif (r1>=interval(1))&(r1 for j=2:iteration-2
if (r1>=interval(j-1))&(r1 choose(1)=j;
break;
end
end
end
%**************************************************************************
%**************************************************************************
%%% 選出第二個(gè)和新點(diǎn)相連接的頂點(diǎn)，注意這兩個(gè)點(diǎn)是不同的，目的是避免重復(fù)邊的出現(xiàn)
r2=rand(1);
if r2 choose(2)=1;
elseif r2>=interval(iteration-2)
choose(2)=iteration-1;
elseif (r2>=interval(1))&(r2 for j=2:iteration-2
if (r2>=interval(j-1))&(r2 choose(2)=j; 
break;
end
end
end

while choose(2)==choose(1)
r2=rand(1);
if r2 choose(2)=1;
elseif r2>=interval(iteration-2)
choose(2)=iteration-1;
elseif (r2>=interval(1))&(r2 for j=2:iteration-2
if (r2>=interval(j-1))&(r2 choose(2)=j; 
break;
end
end
end
end
%**************************************************************************
%**************************************************************************
%%% 選出第三個(gè)和新點(diǎn)相連接的頂點(diǎn)，注意這三個(gè)點(diǎn)是不同的，目的是避免重復(fù)邊的出現(xiàn)
r3=rand(1);
if r3 choose(3)=1;
elseif r3>=interval(iteration-2)
choose(3)=iteration-1;
elseif (r3>=interval(1))&(r3 for j=2:iteration-2
if (r3>=interval(j-1))&(r3 choose(3)=j; 
break;
end
end
end

while (choose(3)==choose(1))|(choose(3)==choose(2))
r3=rand(1);
if r3 choose(3)=1;
elseif r3>=interval(iteration-2)
choose(3)=iteration-1;
elseif (r3>=interval(1))&(r3 for j=2:iteration-2
if (r3>=interval(j-1))&(r3 choose(3)=j;
break;
end
end
end
end
%**************************************************************************
%**************************************************************************
%%% 把新點(diǎn)加入網(wǎng)絡(luò)后，對(duì)鄰接矩陣進(jìn)行相應(yīng)的改變！
%**************************************************************************
%%% 這是在一次循環(huán)下生成的新的鄰接矩陣，下一次循環(huán)就是在這個(gè)鄰接矩陣的基礎(chǔ)上進(jìn)行的！
for k=1:m
adjacent_matrix(iteration,choose(k))=1;
adjacent_matrix(choose(k),iteration)=1;
end
% node_degree=sparse(1,N); % node_degree表示各個(gè)節(jié)點(diǎn)的度數(shù)
for p=1:iteration
%last_element=sparse(iteration,1);
%last_element=cumsum(adjacent_matrix(1:iteration,p));
%node_degree(p)=last_element(iteration);
node_degree(p)=sum(adjacent_matrix(1:iteration,p)); % 這個(gè)循環(huán)的目的是重新給各個(gè)節(jié)點(diǎn)的度賦值
end 
% element_cumsum=sparse(cumsum(adjacent_matrix));
% node_degree=element_cumsum(N,:); 
end 
% 一次最外層循環(huán)的結(jié)束
%**************************************************************************
%**************************************************************************
%**************************************************************************
% element_cumsum=sparse(cumsum(adjacent_matrix)); % element_cumsum的最后一行給出各個(gè)節(jié)點(diǎn)的度數(shù)
% node_degree=element_cumsum(N,:);
number_of_nodes_with_equal_degree=zeros(1,N); % 存儲(chǔ)度相同的頂點(diǎn)的個(gè)數(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é)點(diǎn)的個(gè)數(shù)
% node_degree=element_cumsum(N,:);
end
a_realization_of_distribution=number_of_nodes_with_equal_degree;
for i=1:N
realization_of_distribution(J,i)=a_realization_of_distribution(i); 
end
%%% 循環(huán)完畢之后，清空內(nèi)存，只保留realization_of_distribution的相關(guān)信息，這是唯一有用的數(shù)據(jù)，
%%% 下面的討論僅僅與這個(gè)數(shù)據(jù)有關(guān)
%clear number_of_nodes_with_equal_degree;
%clear element_cumsum;
%clear node_degree;
%clear adjacent_matrix;
% clear 
% clear 
% clear 
% 開始第二次最外層的循環(huán)
%************************************************************************** 
end % 外層循環(huán)的中止
%**************************************************************************
%**************************************************************************
%**************************************************************************
aaa=cumsum(realization_of_distribution);
bb1=aaa(I,:); %%% 譬如，度為3的節(jié)點(diǎn)的個(gè)數(shù)，由于度數(shù)為1,2的節(jié)點(diǎn)的個(gè)數(shù)為0,故可以從度數(shù)為3的節(jié)點(diǎn)個(gè)數(shù)開始計(jì)算
bb2=bb1(m:N);
bbb=bb2/(I*N); %%% 譬如，度為3的節(jié)點(diǎn)的個(gè)數(shù)在網(wǎng)絡(luò)中的比例
K=m:N; %%%% 這是
loglog(K,bbb,'*') % 注意，作圖的時(shí)候，一定要做散點(diǎn)圖
axis([1 N 0.0000001 0.9])
hold on;
y1=2*m^2*K.^(-3);
loglog(K,y1,'r'); % 與平均場(chǎng)結(jié)果進(jìn)行比較 p(k)=2*m^2*k^(-3)

%************************************************************************** 
%%% 第四步：：全部工作結(jié)束
toc; %%% 計(jì)算程序運(yùn)行需要的時(shí)間