%      RandomGraphDemo.M
%   
%  Figures to demonstrate the properties of random graphs
%  
% J. Lunze
% 
% 30.05.2018
% Version v. 6.1.2019
% for 2nd edition: 5.4.2021
%
echo off
clear
close all

%%   Connectivity of random graphs
%
%  determine the lower bound of the connection probability p
%  such that the undirected random graph is probable to be connected
%
figure(1)
subplot(4,1,[1,2])
N1=20;
kmax=150;
Prob=[0.02:0.005:0.4 0.45:0.05:1];
ConnectSum=[];    % average number of connected vertex sets 
                  % in dependence upon the probability p
kp=0;
for p=Prob
    Connect=0;
    %  Determine for the given p the number of connected 
    % subgraphs of kmax realisations
    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
        %  Generate the Laplacian matrix
        D=diag(sum(A));      % degree matrix
        L=D-A;               % Laplacian matrix
        %  determine the number of vanishing eigenvalues
        Eig=eig(L);
        EigSorted=sort(Eig,1,'ascend');
        kzeros=0;
        epsilon=0.001;
        while kzeros<N1 && abs(EigSorted(kzeros+1))<= epsilon
            kzeros=kzeros+1;
        end
        Connect=Connect+kzeros;
    end
    kp=kp+1;
    ConnectSum(kp)=Connect/kmax;
end
plot(Prob, ConnectSum, 'b');
hold on
plot([log(N1)/N1, log(N1)/N1],[0 7],'b:');
latexylabel('$$\bar{s}\;\;$$');
rotateY;
% latexxlabel('$$p$$');
xleer;
axis([0 1 0 7]);
yticks([0  2  4 6]);
%xticks([0 0.2 0.4 0.6 0.8 1]);
%latextext(0.8, 3.4, ['$$N=\;$$', num2str(N1)]); 
latextitle('$${\tt RandomGraphDemo:}$$ Connectivity of random graphs');
hold on
%
%  Determine the probability of the graph to be connected
%
figure(1)
subplot(4,1,[3,4])
Prob=[0.02:0.01:0.4 0.45:0.05:1];
ConnectSum=[];      % probability of the graph to be connected
                    % in dependence upon the probability p
kp=0;
for p=Prob
    Connect=0;
    %  Determine for the given p the number of connected 
    % subgraphs of kmax realisations
    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
        %  Generate the Laplacian matrix
        D=diag(sum(A));      % degree matrix
        L=D-A;               % Laplacian matrix
        %  determine the number of vanishing eigenvalues
        Eig=eig(L);
        EigSorted=sort(Eig,1,'ascend');
        kzeros=0;
        epsilon=0.001;
        while kzeros<N1 && abs(EigSorted(kzeros+1))<= epsilon
            kzeros=kzeros+1;
        end
        if kzeros==1
           Connect=Connect+1;
        end
    end
    kp=kp+1;
    ConnectSum(kp)=Connect/kmax;      % Probability to be connected
end
plot(Prob, ConnectSum, 'b');
%semilogx(Prob, ConnectSum);
hold on
plot([log(N1)/N1, log(N1)/N1],[-0.2 1.1],'b:');
latexylabel('$${\rm Prob}(S=1)$$');
rotateY;
latexxlabel('$$p$$');
axis([0 1 -0.1 1.1]);
yticks([0 0.5 1]);
xticks([0 0.2 0.4 0.6 0.8 1]);
hold on
%
%   the same for N=40
figure(1)
subplot(4,1,[1,2])
N1=45;
ConnectSum=[];
kp=0;
for p=Prob
    Connect=0;
    %  Determine for the given p the number of connected 
    % subgraphs of kmax realisations
    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
        %  Generate the Laplacian matrix
        D=diag(sum(A));      % degree matrix
        L=D-A;               % Laplacian matrix
        %  determine the number of vanishing eigenvalues
        Eig=eig(L);
        EigSorted=sort(Eig,1,'ascend');
        kzeros=0;
        epsilon=0.001;
        while kzeros<N1 && abs(EigSorted(kzeros+1))<= epsilon
            kzeros=kzeros+1;
        end
        Connect=Connect+kzeros;
    end
    kp=kp+1;
    ConnectSum(kp)=Connect/kmax;
end
plot(Prob, ConnectSum, 'b--');
plot([log(N1)/N1, log(N1)/N1],[0 7],'b:');
%latextext(0.9, 2.9, num2str(N1)); 
%latextext(0.8, 2.2, '$$\frac{{\rm ln}(N)}{N}\;$$...'); 
hold off
%
%  Determine the probability of the graph to be connected
%
figure(1)
subplot(4,1,[3,4])
ConnectSum=[];
kp=0;
for p=Prob
    Connect=0;
    %  Determine for the given p the number of connected 
    % subgraphs of kmax realisations
    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
        %  Generate the Laplacian matrix
        D=diag(sum(A));      % degree matrix
        L=D-A;               % Laplacian matrix
        %  determine the number of vanishing eigenvalues
        Eig=eig(L);
        EigSorted=sort(Eig,1,'ascend');
        kzeros=0;
        epsilon=0.001;       % numerical bound for eigenvalue test
        while kzeros<N1 && abs(EigSorted(kzeros+1))<= epsilon
            kzeros=kzeros+1;
        end
        if kzeros==1
           Connect=Connect+1;
        end
    end
    kp=kp+1;
    ConnectSum(kp)=Connect/kmax; % relative frequency to be connected
end
plot(Prob, ConnectSum, 'b--');
plot([log(N1)/N1, log(N1)/N1],[-0.2 1.1],'b:');
hold off

%%  Degree distribution of random graphs
%   Theory: Poisson distribution
figure(2)
subplot(3,1,[1,2]);
N1=51;
p=0.2;
%  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');
hold on
%   compare with realisations of the random graph
%RandomGraphDemo1=rng;                         % save s to restore the generator setting
%save('RandomGraphDemo1', 'RandomGraphDemo1');
load('RandomGraphDemo1');       % Seed (initialisation) of the number generator
rng(RandomGraphDemo1);

%  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 degree distribution
D=sum(A);
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)=DD(degree+1)+1;
            k2=k2+1;
        end
end
DD=DD/N1;
plot(0:N1, DD, 'bo');  
axis([0 N1-1 0 0.2]);
xticks([0 10 20 30 40 50]);
yticks([0 0.1 0.2]);
latexxlabel('$$d$$');
latexylabel('$${\rm Prob}(D_i=d)$$');
rotateY;
%latextext(0.75*N1, 0.18, ['$$N=\;$', num2str(N1)]);
%latextext(0.75*N1, 0.16, ['$$p=\;$', num2str(p)]);
%  determine the average degree of the realisation
dbar=sum(D)/N1;
%plot([dbar,dbar], [0,0.2],'b-.');
plot([p*(N1-1),p*(N1-1)], [0,0.2],'b:');
%latextext(0.75*N1, 0.14, ['$$pN\;$$ -.-']);
%latextext(0.75*N1, 0.12, ['$$\bar{d}\;$$ ...']);
plot([5, 15], [-5 -5], 'k');
hold off

%%  Characteristic path length of random graphs
%   in dependence upon the probability p
figure(3)
subplot(3,1,[1,2]);
N1=40;  % 40
kmax=600;
Prob=[0.07:0.01:0.3 0.35:0.05:1];
Pathlength=[];            % characteristic path length
                    % in dependence upon the probability p
kp=0;
%RandomGraphDemo2=rng;                         % save s to restore the generator setting
%save('RandomGraphDemo2', 'RandomGraphDemo2');
load('RandomGraphDemo2');       % Seed (initialisation) of the number generator
rng(RandomGraphDemo2);
%
for p=Prob
    lbar=0;
    number=0;
    %  Determine for the given p the number of connected 
    % subgraphs of kmax realisations
    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~=0     % if graph is connected
            lbar=lbar+l1;
            number=number+1;   % number of connected graphs
        end
     end
    kp=kp+1;
    Pathlength(kp)=lbar/number;      % Probability to be connected
end
plot(Prob, Pathlength, 'b');
%semilogx(Prob, Pathlength);
hold on
plot([0 1],[1 1],'k:');
plot([log(N1)/N1, log(N1)/N1],[0 7],'b:');
%   Approximation by log N/log(p (N-1))
Lapprox=[];
kp=0;
for p=Prob
    kp=kp+1;
    Lapprox(kp)=2*log10(N1)/log10(p*(N1-1));
end
%plot(Prob, Lapprox, 'b--');   %  approximation
latexylabel('$$\bar{l}$$');
rotateY;
latexxlabel('$$p$$');
axis([0 1 0 7]);
yticks([ 1 2 3 5 7]);
xticks([0  0.4 0.6 0.8 1]);
latextitle('$${\tt RandomGraphDemo:}$$ Characteristic path length of random graphs');
latextext(0.8, 3.5, ['$$N=\;$$', num2str(N1)]);
plot(0.22, 2, 'bo');
plot([0 0.22],[2 2],'b:', [0.22 0.22], [-1 2], 'b:');
latextext(0.15,-0.6,'$$\bar{p}=0.22$$');
hold off


%%   Random graph - illustration of characteristic properties
%
N1=15;
%RandomGraphDemo3=rng;                         % save s to restore the generator setting
%save('RandomGraphDemo3', 'RandomGraphDemo3');
load('RandomGraphDemo3');       % Seed (initialisation) of the number generator
rng(RandomGraphDemo3);

%
figure(4)
p=0.05;
p1=p/2;
%  Generation of the adjacency matrix of a random graph
A=ceil(rand(N1,N1)-(1-p1)*ones(N1,N1)); 
%  make the matrix symmetric for undirected graphs
A=floor(A+A');
% set the main diagonal elements zero
A=A-diag(diag(A)); 
% sort the connected components by size
A=findConnComp(A);
G=graph(A);
plot(G, 'Layout', 'circle', 'LineWidth', 2, 'MarkerSize', 8);  
yticks([]);
xticks([]);
axis('square');
latextext(0.8, 1.3, num2str(p));
%
figure(5)
p=0.1;
p1=p/2;
A=ceil(rand(N1,N1)-(1-p1)*ones(N1,N1)); 
A=floor(A+A');
A=A-diag(diag(A)); 
A=findConnComp(A);
G=graph(A);
plot(G, 'Layout', 'circle', 'LineWidth', 2, 'MarkerSize', 8);  
yticks([]);
xticks([]);
axis('square');
latextext(0.8, 1.3, num2str(p));
%
figure(6)
p=0.15;
p1=p/2;
A=ceil(rand(N1,N1)-(1-p1)*ones(N1,N1)); 
A=floor(A+A');
A=A-diag(diag(A)); 
A=findConnComp(A);
G=graph(A);
plot(G, 'Layout', 'circle', 'LineWidth', 2, 'MarkerSize', 8);  
yticks([]);
xticks([]);
axis('square');
latextext(0.8, 1.3, num2str(p));
%
figure(7)
p=0.2;
p1=p/2;
A=ceil(rand(N1,N1)-(1-p1)*ones(N1,N1)); 
A=floor(A+A');
A=A-diag(diag(A)); 
A=findConnComp(A);
G=graph(A);
plot(G, 'Layout', 'circle', 'LineWidth', 2, 'MarkerSize', 8);  
yticks([]);
xticks([]);
axis('square');
latextext(0.8, 1.3, num2str(p));
%
figure(8)
p=0.3;
p1=p/2;
A=ceil(rand(N1,N1)-(1-p1)*ones(N1,N1)); 
A=floor(A+A');
A=A-diag(diag(A)); 
A=findConnComp(A);
G=graph(A);
plot(G, 'Layout', 'circle', 'LineWidth', 2, 'MarkerSize', 8);  
yticks([]);
xticks([]);
axis('square');
latextext(0.8, 1.3, num2str(p));
%
figure(9)
p=0.4;
p1=p/2;
A=ceil(rand(N1,N1)-(1-p1)*ones(N1,N1)); 
A=floor(A+A');
A=A-diag(diag(A)); 
A=findConnComp(A);
G=graph(A);
plot(G, 'Layout', 'circle', 'LineWidth', 2, 'MarkerSize', 8);  
yticks([]);
xticks([]);
axis('square');
latextext(0.8, 1.3, num2str(p));
%


%%  Degree distribution of the example graph
%   Compare the Poisson distribution with the distribution of the example
%   graph
figure(10)
subplot(3,1,[1,2]);
N1=12;
lambda=3.4;
%  Poisson distribution
Prob=[];
for d=0:N1
%    Prob(d+1)=nchoosek(N1,d)*lambda^d*(1-lambda)^(N1-d);     % nchoosek determines the binomial coefficient
     Prob(d+1)=exp(-lambda)*lambda^d/factorial(d);
end
%   draw the bar chart of the example graph
H=[0 3 6 2 0 1]/N1;   % relative frequency of vertex degree
D=[1 2 3 4 5 6];       % vertex degree
for k1=1:length(D)
   fill([D(k1)-0.3 D(k1)-0.3 D(k1)+0.3 D(k1)+0.3 D(k1)-0.3],...
        [0 H(k1), H(k1), 0, 0], [122, 187, 239]/250); 
    hold on
end
%   draw the Poisson distribution 
plot(0:N1, Prob, 'bo',...
            'MarkerEdgeColor','b',...
            'MarkerFaceColor','b');
plot(0:N1, Prob, 'b--');
latexylabel('$${\rm Prob}(D=d)$$');
rotateY;
latexxlabel('$$d$$');
axis([0 N1-1 0 0.6]);
yticks([ 0 0.5 1]);
xticks([1 2 3 4 5 6 7 8 9 10 11 12]);
latextitle('$${\tt RandomGraphDemo:}$$ Degree distribution');
hold off



%%  Characteristic path length of random graphs
%   in dependence upon the vertex number N1
%
figure(11)
subplot(3,1,[1,2]);
dbar=15;      % mean vertex degree  dbar = p(N-1)
kmax=25;
Nmax=1000;
Vertexnumber=50:50:Nmax;
Pathlength=zeros(size(Vertexnumber));      % characteristic path length
                         % in dependence upon the number N1 of vertices
kp=0;
%
for N1=Vertexnumber
    p=dbar/(N1-1);
    lbar=0;
    number=0;
    %  Determine for the given p, N1 the characteristic path length
    %  of kmax realisations
    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~=0     % if graph is connected
            lbar=lbar+l1;
            number=number+1;   % number of connected graphs
        end
     end
    kp=kp+1;
    Pathlength(kp)=lbar/number;    % Mean of the char. path lengths
end
plot(Vertexnumber, Pathlength, 'b');
hold on
%   Approximation by log N/log(p (N-1))
Lapprox=zeros(size(Pathlength));
kp=0;
for N1=Vertexnumber
    kp=kp+1;
    Lapprox(kp)=log10(N1)/log10(dbar);
end
plot(Vertexnumber, Lapprox, 'b--');
latexylabel('$$l_{\rm max}$$');
rotateY;
latexxlabel('$$N$$');
axis([0 Nmax 0 4]);
yticks([ 1 2 3 4]);
xticks([200 400 600 800 1000]);
latextitle('$${\tt RandomGraphDemo:}$$ Diameter of random graphs');
latextext(800, 3.5, ['$$\bar{d}=\;$$', num2str(dbar)]);
hold off


%%
epsfigc('RandomGraphDemo');



