%      ScaleFreeDemo.M
%   
%  Figures to demonstrate the properties of scale-free networks
%  
% J. Lunze
% 
% 08.06.2018
% Version v. 6.1.2019
%
echo off
clear
close all

%%   Comparison of exponential distribution with powerlaw distribution
%
figure(1)
subplot(3,1,[1 2])
dmax=20;       %  maximum degree considered
D=1:dmax;
%     exponential probability distribution
lambda=1.7;
Exponential=exp(-lambda*D);
Exponential=Exponential/sum(Exponential);  % get sum of probabilities = 1
plot(D, Exponential, 'bo');
hold on
plot(D, Exponential, 'b--');
%      powerlaw distribution
alpha=1.2;
Powerlaw=D.^(-alpha);  % determine elementwise power
Powerlaw=Powerlaw/sum(Powerlaw);     % get sum of probabilities = 1
plot(D, Powerlaw, 'bo',...
            'MarkerEdgeColor','b',...
            'MarkerFaceColor','b');
hold on
plot(D, Powerlaw, 'b');
%
latextitle('$${\tt ScaleFreeDemo:}$$ Exponential vs. powerlaw distribution');
latexxlabel('$$d$$');
latexylabel('$${\rm Prob}(D_i=d)$$');
rotateY;
xticks([1 5 10 15 20]);
yticks([0 0.2 0.4]);
plot([0 dmax],[0 0],'k:');
axis([0 dmax -0.02 0.4]);
hold off
%   probability to have degree d>5
ProbExp=0;
for k1=6:dmax
    ProbExp=ProbExp+Exponential(k1);
end
ProbExp
ProbPower=0;
for k1=6:dmax
    ProbPower=ProbPower+Powerlaw(k1);
end
ProbPower

%%  Properties of scale-free networks
%
figure(2)
Nmax=20;   % number of vertices in the final graph
N0=6;      % initial vertex number
d0=3;    % degree of new vertices
fig1=N0;    % graph to be shown in a figure: fig1 = number of vertices connected
fig2=N0+8;
fig3=Nmax;
% dbar=3;    % degree of new vertices
% A0=[0 1 0 1 1 0;   % adjacency matrix of initial graph
%     1 0 1 0 0 1;
%     0 1 0 1 1 0;
%     1 0 1 0 0 1;
%     1 0 1 0 0 1;
%     0 1 0 1 1 0];
% Vertices=[1 2 3 4 5 6 1 2 3 4 5 6 1 2 3 4 5 6]; % vector with vertices
%                            % with multiplicity according to their degree
A0=[0 1 0 0 1 0;   % adjacency matrix of initial graph
    1 0 0 0 0 1;
    0 0 0 0 1 0;
    0 0 0 0 0 1;
    1 0 1 0 0 1;
    0 1 0 1 1 0];
Vertices=[1 2 3 4 5 6 1 2 3 4 5 6 1 2 3 4 5 6]; % vector with vertices
                           % with multiplicity according to their degree
mk=sum(sum(A0));     % 2 * number of edges in the current graph
                     %  = length(Vertices)
A=zeros(Nmax, Nmax);
A(1:N0, 1:N0)=A0;
   G=graph(A);
   plot(G, 'Layout', 'circle', 'LineWidth', 2, 'MarkerSize', 8);  
   yticks([]);
   xticks([]);
   axis('square');
%
for i1=N0+1:Nmax
    Prob=rand([1,mk]);
    [Prob,Index]=sort(Prob,'descend');  % begins with highest probability
    Endpoints(1)=Vertices(Index(1));   % first new endpoint of edges
    k1=1;    % number of vertices found
    k2=1;
    while k1<d0
        k2=k2+1;
        if ~ismember(Vertices(k2),Endpoints)
            k1=k1+1;
            Endpoints(k1)=Vertices(k2);
        end
    end
   for k1=1:d0
       A(i1,Endpoints(k1))=1;   % Vertices(Index(k1)) = vertex number
                                      % that belongs to the Index(k1)
       A(Endpoints(k1),i1)=1;
       Vertices=[Vertices,Endpoints(k1),i1];   % add the two vertices
   end
   mk=mk+2*d0;
   if i1==fig2
       figure(3)
       G=graph(A);
       plot(G, 'Layout', 'circle', 'LineWidth', 2, 'MarkerSize', 8);  
       yticks([]);
       xticks([]);
       axis('square');
   end
   if i1==fig3
       figure(4)
       G=graph(A);
       plot(G, 'Layout', 'circle', 'LineWidth', 2, 'MarkerSize', 8);  
       yticks([]);
       xticks([]);
       axis('square');
   end
end



%%  Degree distribution
%
figure(5)
subplot(3,1,[1,2]);
Nmax=150;   % number of vertices in the final graph
N0=6;      % initial vertex number
d0=3;    % degree of new vertices
%ScaleFreeDemo1=rng;                         % save s to restore the generator setting
%save('ScaleFreeDemo1', 'ScaleFreeDemo1');
load('ScaleFreeDemo1');       % Seed (initialisation) of the number generator
rng(ScaleFreeDemo1);

% A0=[0 1 0 1 1 0;   % adjacency matrix of initial graph
%     1 0 1 0 0 1;
%     0 1 0 1 1 0;
%     1 0 1 0 0 1;
%     1 0 1 0 0 1;
%     0 1 0 1 1 0];
% Vertices=[1 2 3 4 5 6 1 2 3 4 5 6 1 2 3 4 5 6]; % vector with vertices
%                            % with multiplicity according to their degree
A0=[0 1 0 0 1 0;   % adjacency matrix of initial graph
    1 0 0 0 0 1;
    0 0 0 0 1 0;
    0 0 0 0 0 1;
    1 0 1 0 0 1;
    0 1 0 1 1 0];
Vertices=[1 2 3 4 5 6 1 2 3 4 5 6 1 2 3 4 5 6]; % vector with vertices
                           % with multiplicity according to their degree
mk=sum(sum(A0));     % 2 * number of edges in the current graph
                     %  = length(Vertices)
A=zeros(Nmax, Nmax);
A(1:N0, 1:N0)=A0;
for i1=N0+1:Nmax
    Prob=rand([1,mk]);
    [Prob,Index]=sort(Prob,'descend');  % begins with highest probability
    Endpoints(1)=Vertices(Index(1));   % first new endpoint of edges
    k1=1;    % number of vertices found
    k2=1;
    while k1<d0
        k2=k2+1;
        if ~ismember(Vertices(k2),Endpoints)
            k1=k1+1;
            Endpoints(k1)=Vertices(k2);
        end
    end
   for k1=1:d0
       A(i1,Endpoints(k1))=1;   % Vertices(Index(k1)) = vertex number
                                      % that belongs to the Index(k1)
       A(Endpoints(k1),i1)=1;
       Vertices=[Vertices,Endpoints(k1),i1];   % add the two vertices
   end
   mk=mk+2*d0;
end
%  Evaluation of the network: bar chart of vertex degrees
%  determine the degree distribution (from RandomGraphDemo.m)
D=sum(A);
D=sort(D, 'ascend');
DD=zeros(Nmax+1,1);      % degree distribution: 
                       % DD(k1) = number of vertices with degree k1-1
k2=1;
for degree=0:Nmax
        while k2 <= Nmax && D(k2)==degree 
            DD(degree+1)=DD(degree+1)+1;
            k2=k2+1;
        end
end
DD=DD/Nmax;
%  plot the non-zero values
for k1=1:Nmax+1
    if DD(k1)~=0
        plot(k1-1, DD(k1), 'bo');
        hold on
    end
end
plot([0 Nmax],[0 0],'k:');
plot(0:Nmax, DD, 'b:');  
axis([0 90 -0.03 0.6]);
xticks([0 10 30 50 70 90]);
yticks([0 0.2 0.4 0.6]);
latexxlabel('$$d$$');
latexylabel('$${\rm Prob}(D_i=d)$$');
rotateY;
%  determine the average degree of the realisation
d0=sum(D)/Nmax;
%plot([p*Nmax,p*Nmax], [0,0.2],'b:');
latextitle('$${\tt ScaleFreeDemo:}$$ Degree distribution');
latextext(70, 0.5, ['$$N=\;$', num2str(Nmax)]);
hold off


%%  smoothing probability distribution
%  take the "inverse" rule: connect to vertices with low degree
%
figure(6)
subplot(3,1,[1,2]);
Nmax=150;   % number of vertices in the final graph
N0=6;      % initial vertex number
d0=3;    % degree of new vertices
A0=[0 1 0 0 1 0;   % adjacency matrix of initial graph
    1 0 0 0 0 1;
    0 0 0 0 1 0;
    0 0 0 0 0 1;
    1 0 1 0 0 1;
    0 1 0 1 1 0];
Vertices=[1 2 3 4 5 6 1 2 3 4 5 6 1 2 3 4 5 6]; % vector with vertices
                           % with multiplicity according to their degree
mk=sum(sum(A0));     % 2 * number of edges in the current graph
                     %  = length(Vertices)
A=zeros(Nmax, Nmax);
A(1:N0, 1:N0)=A0;
%
for i1=N0+1:Nmax
    Prob=rand([1,mk]);
    [Prob,Index]=sort(Prob,'ascend');  % begins with smallest probability
    % Does this ordering ensure that the "inverse" prob distribution is
    % applied?
    Endpoints(1)=Vertices(Index(1));   % first new endpoint of edges
    k1=1;    % number of vertices found
    k2=1;
    while k1<d0
        k2=k2+1;
        if ~ismember(Vertices(k2),Endpoints)
            k1=k1+1;
            Endpoints(k1)=Vertices(k2);
        end
    end
   for k1=1:d0
       A(i1,Endpoints(k1))=1;   % Vertices(Index(k1)) = vertex number
                                      % that belongs to the Index(k1)
       A(Endpoints(k1),i1)=1;
       Vertices=[Vertices,Endpoints(k1),i1];   % add the two vertices
   end
   mk=mk+2*d0;
end
%  Evaluation of the network: bar chart of vertex degrees
%  determine the degree distribution (from RandomGraphDemo.m)
D=sum(A);
D=sort(D, 'ascend');
DD=zeros(Nmax+1,1);      % degree distribution: 
                       % DD(k1) = number of vertices with degree k1-1
k2=1;
for degree=0:Nmax
        while k2 <= Nmax && D(k2)==degree 
            DD(degree+1)=DD(degree+1)+1;
            k2=k2+1;
        end
end
DD=DD/Nmax;
%  plot the non-zero values
for k1=1:Nmax+1
    if DD(k1)~=0
        plot(k1-1, DD(k1), 'bo');
        hold on
    end
end
plot([0 Nmax],[0 0],'k:');
plot(0:Nmax, DD, 'b:');  
axis([0 90 -0.03 1.1*max(DD)]);
xticks([0 10 30 50 70 90]);
yticks([0 0.2 0.4 0.6]);
latexxlabel('$$d$$');
latexylabel('$${\rm Prob}(D_i=d)$$');
rotateY;
%  determine the average degree of the realisation
d0=sum(D)/Nmax;
%plot([p*Nmax,p*Nmax], [0,0.2],'b:');
latextitle('$${\tt ScaleFreeDemo:}$$ Degree distribution for the inverse rule');
latextext(70, 0.5, ['$$N=\;$', num2str(Nmax)]);
hold off


%%  Characteristic path length
%
figure(7)
subplot(3,1,[1,2]);
Nmax=90;   % number of vertices in the final graph
N0=6;      % initial vertex number
d0=3;    % degree of new vertices
kmax=10;
A0=[0 1 0 0 1 0;   % adjacency matrix of initial graph
    1 0 0 0 0 1;
    0 0 0 0 1 0;
    0 0 0 0 0 1;
    1 0 1 0 0 1;
    0 1 0 1 1 0];
Vertices=[1 2 3 4 5 6 1 2 3 4 5 6 1 2 3 4 5 6]; % vector with vertices
                           % with multiplicity according to their degree
mk=sum(sum(A0));     % 2 * number of edges in the current graph
                     %  = length(Vertices)
kp=0;
PathLength=zeros(Nmax-2*N0+1,1);
for N1=2*N0:3:Nmax
    kp=kp+1;
    Length=0;
    for k2=1:kmax
        %  generate a scale-free network
        A=zeros(Nmax, Nmax);
        A(1:N0, 1:N0)=A0;
        for i1=N0+1:Nmax
            Prob=rand([1,mk]);
            [Prob,Index]=sort(Prob,'descend');  % begins with highest probability
            Endpoints(1)=Vertices(Index(1));   % first new endpoint of edges
            k1=1;    % number of vertices found
            k2=1;
            while k1<d0
                k2=k2+1;
                if ~ismember(Vertices(k2),Endpoints)
                    k1=k1+1;
                    Endpoints(k1)=Vertices(k2);
                end
            end
           for k1=1:d0
               A(i1,Endpoints(k1))=1;   % Vertices(Index(k1)) = vertex number
                                              % that belongs to the Index(k1)
               A(Endpoints(k1),i1)=1;
               Vertices=[Vertices,Endpoints(k1),i1];   % add the two vertices
           end
           mk=mk+2*d0;
        end
        Length=Length+charPathLength(A);
    end
    PathLength(kp,1)=Length/kmax;
end
plot(2*N0:Nmax, PathLength, 'b');  
axis([0 Nmax 0 3]);
xticks([0 10 20 30 40 50 60]);
yticks([0 1 2 3 ]);
latexxlabel('$$N$$');
latexylabel('$$\bar{l}$$');
rotateY;
latextitle('$${\tt ScaleFreeDemo:}$$ Characteristic path length');
latextext(70, 0.5, ['$$\bar{d}=\;$', num2str(d0)]);
hold off





%%
epsfigc('ScaleFreeDemo');