%      SmallWorldDemo.M
%   
%  Figures to demonstrate the properties of small-world networks
%  
% J. Lunze
% 
% 02.06.2018
% Version v. 2.7.2018
% for 2nd edition: 6.4.2021
%
echo off
clear
close all

%%  Small-world model
%   starting point: undirected double ring
%   pr - re-wiring probability
%   lbar - characteristic path length
figure(1)
subplot(3,1,[1,2]);
N1=200;      % 60
kNeigh=4;   % number of neighbours in the regular graph
kmax=25;    % 250
Prob=[0.001 0.002 0.005 0.01:0.01:0.2, 0.22:0.1:1 1];  % re-connection probability
%  adjacency matrix of the double ring
ARing=diag(ones(N1-1,1),-1)+diag(ones(N1-1,1),1);
ARing=ARing+diag(ones(N1-2,1),-2)+diag(ones(N1-2,1),2);
ARing(1,N1)=1;
ARing(N1,1)=1;
ARing(2,N1)=1;
ARing(N1,2)=1;
ARing(N1-1,1)=1;
ARing(1,N1-1)=1;
Noe=sum(sum(ARing))/2;  % number of undirected edges
%
Pathlength=[];
Clustercoeff=[];
ClusterApprox=[];
kp=1;
Pathlength(1)=charPathLength(ARing);  % regular graph
Clustercoeff(1)=clustercoeff(ARing);
ClusterApprox(1)=3*(kNeigh-1)/2/(2*kNeigh-1);
for pr=Prob
    pr
    lbar=0;
    knew=0;
    cluster=0;
    for k3=1:kmax
        Rewire=ceil(rand(Noe,1)-(1-pr)*ones(Noe,1)); % vector of elements in {0,1} 
                           % determining which edge should be re-connected
        NewDest=ceil(N1*rand(Noe,1)); % vector of new destinations of the edge
        %   re-connect the edges
        %   kedge - edge number
        kedge=0;
        ARingNew=zeros(N1,N1);
        for k1=1:N1
            for k2=k1+1:N1
                if ARing(k1,k2)==1
                    kedge=kedge+1;
                    if Rewire(kedge,1)==1 && ARing(NewDest(kedge,1),k2)~=1 ...
                            && NewDest(kedge,1)~=k2
                       ARingNew(NewDest(kedge,1),k2)=1;
                       ARingNew(k2,NewDest(kedge,1))=1;
                    else
                       ARingNew(k1,k2)=1;
                       ARingNew(k2,k1)=1;
                    end
                end
            end
        end
        %   determine the characteristic path length
        lnew=charPathLength(ARingNew);
        if lnew ~=N1    % otherwise the graph is not connected
            knew=knew+1;
            lbar=lbar+lnew;
        end
        cluster=cluster+clustercoeff(ARingNew);
    end
    kp=kp+1;
    Pathlength(kp)=lbar/knew;
    Clustercoeff(kp)=cluster/kmax;
    ClusterApprox(kp)=3*(kNeigh-1)/2/(2*kNeigh-1)*(1-pr)^3;
end
semilogx([0 1],[0 0],'k');    %   to create a figure with logarithmic axis
hold on
fill([0.01 0.08 0.08 0.01 0.01],[0.01 0.01 0.99 0.99 0.01],[0.94 0.94 0.94],'EdgeColor','none');
semilogx([0 Prob], Pathlength/Pathlength(1), 'b');
semilogx([0 Prob], Clustercoeff/max(Clustercoeff), 'b');
semilogx([0 Prob], ClusterApprox/max(ClusterApprox), 'b--');
% plot([0 Prob], Pathlength/max(Pathlength), 'b');
% hold on
% plot([0 Prob], Clustercoeff/max(Clustercoeff), 'b');
latextitle('$${\tt SmallWorldDemo:}$$ Path length of small-world network');
latexxlabel('$$p_{\rm r}$$');
latexylabel('$$\frac{\bar{l}}{\bar{l}(0)}, \frac{c}{c(0)}$$');
rotateY;
xticks([0.001 0.01 0.1 1]);
yticks([0 0.5 1]);
axis([0.001 1 0 1]);     
latextext(0.2, 0.85, ['$$N=\;$$', num2str(N1)]);
hold off
%
%  Determination of the characteristic path length for the corresponding
%  random graph
%  from RandomGraphDemo.m
mmax=(N1-1)*N1/2;   % maximum number of edges
mactual=kNeigh*N1/2; % actual number of edges
p=mactual/mmax;
kmax=30;
%
lbar=0;            % characteristic path length
number=0;
for k1=1:kmax
        %  Generation of the adjacency matrix of a random graph
        A=ceil(rand(N1,N1)-(1-p)*ones(N1,N1)); 
        %  make the matrix symmetric for undirected graphs
        for k3=1:N1
            for k4=k3+1:N1
                A(k3, k4)=A(k4, k3);
            end
        end
        A=A-diag(diag(A));     % set the main diagonal elements zero
        %  determine the characteristic path length
        l1=charPathLength(A);
        if l1~=N1     % if graph is connected
            lbar=lbar+l1;
            number=number+1;   % number of connected graphs
        end
end
lbar=lbar/number      % characteristic path length
 

%%  Small-world network illustration in three figures
%   starting point: undirected double ring
%   pr - re-wiring probability
%   lbar - characteristic path length
figure(2)
N1=17;
kNeigh=4;   % number of neighbours in the regular graph
pr=0;
%  adjacency matrix of the double ring
ARing=diag(ones(N1-1,1),-1)+diag(ones(N1-1,1),1);
ARing=ARing+diag(ones(N1-2,1),-2)+diag(ones(N1-2,1),2);
ARing(1,N1)=1;
ARing(N1,1)=1;
ARing(2,N1)=1;
ARing(N1,2)=1;
ARing(1,N1-1)=1;
ARing(N1-1,1)=1;
G=graph(ARing);
plot(G, 'Layout', 'circle', 'LineWidth', 2, 'MarkerSize', 8);  
yticks([]);
xticks([]);
axis('square');
lbar=charPathLength(ARing);
latextext(0.4, -1.3, ['$$\bar{l}=\;$$' num2str(lbar,3)]);
latextext(0.4, 1.3, ['$$p_{\rm r}=\;$$' num2str(pr)]);
%
figure(3)
pr=0.09;
%SmallWorldDemo1=rng;                         % save s to restore the generator setting
%save('SmallWorldDemo1', 'SmallWorldDemo1');
load('SmallWorldDemo1');       % Seed (initialisation) of the number generator
rng(SmallWorldDemo1);
Noe=sum(sum(ARing))/2;  % number of undirected edges
%
Rewire=ceil(rand(Noe,1)-(1-pr)*ones(Noe,1)); % vector of elements in {0,1} 
%   determine which edge should be re-wired
NewDest=ceil(N1*rand(Noe,1)); % vector of new destinations of the edge
%   re-wire the edges
kedge=0;       %   kedge - edge number
ARingNew=zeros(N1,N1);
for k1=1:N1
    for k2=k1+1:N1
        if ARing(k1,k2)==1
            kedge=kedge+1;
            if Rewire(kedge,1)==1 && ARing(NewDest(kedge,1),k2)~=1 ...
                    && NewDest(kedge,1)~=k2
                ARingNew(NewDest(kedge,1),k2)=1;
                ARingNew(k2,NewDest(kedge,1))=1;
            else
                ARingNew(k1,k2)=1;
                ARingNew(k2,k1)=1;
            end
        end
    end
end
G=graph(ARingNew);
plot(G, 'Layout', 'circle', 'LineWidth', 2, 'MarkerSize', 8);  
yticks([]);
xticks([]);
axis('square');
lbar=charPathLength(ARingNew);
latextext(0.4, -1.3, ['$$\bar{l}=\;$$' num2str(lbar,3)]);
latextext(0.4, 1.3, ['$$p_{\rm r}=\;$$' num2str(pr)]);
%
figure(4)
pr=0.22;
%
Rewire=ceil(rand(Noe,1)-(1-pr)*ones(Noe,1)); % vector of elements in {0,1} 
%   re-wire the edges
kedge=0;       %   kedge - edge number
ARingNew=zeros(N1,N1);
for k1=1:N1
    for k2=k1+1:N1
        if ARing(k1,k2)==1
            kedge=kedge+1;
            if Rewire(kedge,1)==1 && ARing(NewDest(kedge,1),k2)~=1 ...
                    && NewDest(kedge,1)~=k2
                ARingNew(NewDest(kedge,1),k2)=1;
                ARingNew(k2,NewDest(kedge,1))=1;
            else
                ARingNew(k1,k2)=1;
                ARingNew(k2,k1)=1;
            end
        end
    end
end
G=graph(ARingNew);
plot(G, 'Layout', 'circle', 'LineWidth', 2, 'MarkerSize', 8);  
yticks([]);
xticks([]);
axis('square');
lbar=charPathLength(ARingNew);
latextext(0.4, -1.3, ['$$\bar{l}=\;$$' num2str(lbar,3)]);
latextext(0.4, 1.3, ['$$p_{\rm r}=\;$$' num2str(pr)]);


%%  Degree distribution in dependence upon the re-connection probability
kFig=5;     % figures 5 - 8
N1=60;      % 60
kNeigh=4;   % number of neighbours in the regular graph
Prob=[0.01 0.1 0.4 1];  % re-connection probability
%  adjacency matrix of the double ring
ARing=diag(ones(N1-1,1),-1)+diag(ones(N1-1,1),1);
ARing=ARing+diag(ones(N1-2,1),-2)+diag(ones(N1-2,1),2);
ARing(1,N1)=1;
ARing(N1,1)=1;
ARing(2,N1)=1;
ARing(N1,2)=1;
ARing(N1-1,1)=1;
ARing(1,N1-1)=1;
Noe=sum(sum(ARing))/2;  % number of undirected edges
%
%SmallWorldDemo2=rng;                         % save s to restore the generator setting
%save('SmallWorldDemo2', 'SmallWorldDemo2');
load('SmallWorldDemo2');       % Seed (initialisation) of the number generator
rng(SmallWorldDemo2);
for pr=Prob
    Rewire=ceil(rand(Noe,1)-(1-pr)*ones(Noe,1)); % vector of elements in {0,1} 
                      % determining which edge should be re-connected
    NewDest=ceil(N1*rand(Noe,1)); % vector of new destinations of the edge
    %   re-connect the edges
    %   kedge - edge number
    kedge=0;
    ARingNew=zeros(N1,N1);
        for k1=1:N1
            for k2=k1+1:N1
                if ARing(k1,k2)==1
                    kedge=kedge+1;
                    if Rewire(kedge,1)==1 && ARing(NewDest(kedge,1),k2)~=1 ...
                            && NewDest(kedge,1)~=k2
                       ARingNew(NewDest(kedge,1),k2)=1;
                       ARingNew(k2,NewDest(kedge,1))=1;
                    else
                       ARingNew(k1,k2)=1;
                       ARingNew(k2,k1)=1;
                    end
                end
            end
        end
    %  determine the degree distribution
    D=sum(ARingNew);
    D=sort(D, 'ascend');
    DD=zeros(N1+1,1);      % degree distribution: 
                           % DD(k1) = number of vertices with degree k1-1
    k2=1;
    for degree=0:N1
        while k2 <= N1 && D(k2)==degree 
            DD(degree+1,1)=DD(degree+1,1)+1;
            k2=k2+1;
        end
    end
    DD=DD/N1;
    figure(kFig);
    subplot(3,1,[1,2]);
    kFig=kFig+1;
    plot(0:N1, DD, 'bo');  
    hold on
    plot(0:N1, DD, 'b:');
    %  degree distribution of a similar random network
    mmax=(N1-1)*N1/2;      % maximum number of edges in the network
    p=Noe/mmax;            % approximate connection probability
    %  draw the poisson distribution
    Prob=[];
    for k1=0:N1
        Prob(k1+1)=nchoosek(N1,k1)*p^k1*(1-p)^(N1-k1);     
                            % nchoosek determines the binomial coefficient
    end
    plot(0:N1, Prob, 'b--');
    latextext(15, 0.5*max(DD), ['$$p_{\rm r}=\;$$', num2str(pr)]);
    axis([0 20 0 1.1*max(DD)]);
    xticks([0 2 4 6 10 15 20]);
    yticks([0 0.2 0.4 0.6 0.8 1]);
    latexxlabel('$$d$$');
    latextitle('$${\tt SmallWorldDemo:}$$ Degree distribution of small-world network');
    hold off
end
% figure(kFig-1)
% latexylabel('$${\rm Prob}(D_i=d)$$');
%rotateY;


%%  Small-world model
%   starting point: undirected double ring
%   show the reduction to 1/2 pathlength by 4 shortcuts
%   independently of N
figure(9)
subplot(3,1,[1,2]);
Nmax=250;      % 60
N1set=[10:4:Nmax];
kNeigh=4;   % number of neighbours in the regular graph
kmax=20;    % number of realisations considered
%
relPathlength=[];    % relative path length=char. path length/path length of ring graph
kp=0;
for N1=N1set
    %  adjacency matrix of the double ring
    ARing=diag(ones(N1-1,1),-1)+diag(ones(N1-1,1),1);
    ARing=ARing+diag(ones(N1-2,1),-2)+diag(ones(N1-2,1),2);
    ARing(1,N1)=1;
    ARing(N1,1)=1;
    ARing(2,N1)=1;
    ARing(N1,2)=1;
    ARing(N1-1,1)=1;
    ARing(1,N1-1)=1;
    Noe=sum(sum(ARing))/2;  % number of undirected edges
    %
    LengthRing=charPathLength(ARing);  % regular graph
    lbar=0;
    knew=0;
    for k3=1:kmax
        Rewire=(rand(Noe,1)); % vector of probability values in {0,1} 
                           % determining which edge should be re-connected
        NewDest=ceil(N1*rand(Noe,1)); % vector of new destinations of the edge
        %   re-connect the edges
        EdgeReconn=[];    % index of edges to be reconnected
        EdgeDest=[];      % index of new vertex
        for k4=1:4        
            [Max,imax]=max(Rewire); % find index for the edge to be reconnected
            EdgeReconn=[EdgeReconn imax];
            EdgeDest=[EdgeDest NewDest(imax)];
            Rewire(imax)=0;         % reset to find the next probable edge
        end
        %  find edges to be reconnected
        ARingNew=ARing;
        for k4=1:4
            kedge=0;     % enumerate the edges
            for k1=1:N1
                for k2=k1+1:N1
                    if ARingNew(k1,k2)==1
                        kedge=kedge+1;
                        if kedge==EdgeReconn(k4)   % edge found
                           ARingNew(k1,k2)=0;
                           ARingNew(k1,EdgeDest(k4))=1;
                           ARingNew(k2,k1)=0;
                           ARingNew(EdgeDest(k4),k1)=1;
                        end
                    end
                end
            end
        end
        %   determine the characteristic path length
        lnew=charPathLength(ARingNew);
        if lnew ~=N1    % otherwise the graph is not connected
            knew=knew+1;
            lbar=lbar+lnew;
        end
    end
    kp=kp+1;
    relPathlength(kp)=lbar/knew/LengthRing;
end
plot(N1set, relPathlength, 'b');
hold on
plot([0 Nmax],[0.5 0.5],'b--');
latextitle('$${\tt SmallWorldDemo:}$$ Path length of small-world network');
latexxlabel('$$N$$');
latexylabel('$$\frac{\bar{l}}{\bar{l}_{\rm ring}}$$');
rotateY;
xticks([10 50 100 150 200 250]);
yticks([0 0.5 1]);
axis([0 Nmax 0 1.1]);     
hold off

 

%%
epsfigc('SmallWorldDemo');

        
        