%      RandomGraphs.M
%   
%  Figures to demonstrate the properties of random graphs
%  
% J. Lunze
% 
% 1.08.2023
% Version v. 4.8.2023
%
echo off
clear
close all



%%  Distribution of random graphs
%   Theory: Binomial distribution
figure(1)
subplot(3,1,[1,2]);
figure(2)
subplot(3,1,[1,2]);
N1=25;
mmax=(N1-1)*N1/2;
for p=[0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9]
%  draw the Binomial distribution
Prob=[];
Sum=[];
Prob(1)=(1-p)^mmax;
Sum(1)=Prob(1);
for k1=1:mmax
    Prob(k1+1)=nchoosek(mmax,k1)*p^k1*(1-p)^(mmax-k1);     
                     % nchoosek determines the binomial coefficient
    Sum(k1+1)=Sum(k1)+Prob(k1+1);
end
figure(1)
plot(0:mmax, Prob, 'b');
hold on
%   expected value for the number of edges
Em=p*mmax;
plot([Em Em],[0 max(Prob)],'b--');
figure(2)
plot(0:mmax, Sum, 'b');
hold on
end
figure(1)
latextext(5, 0.1, '$$p=$$');
latextext(20, 0.09, '$$0.1$$');
latextext(50, 0.09, '$$0.2$$');
latextext(80, 0.09, '$$0.3$$');
latextext(110, 0.09, '$$0.4$$');
latextext(140, 0.09, '$$0.5$$');
latextext(170, 0.09, '$$0.6$$');
latextext(200, 0.09, '$$0.7$$');
latextext(230, 0.09, '$$0.8$$');
latextext(260, 0.09, '$$0.9$$');
axis([0 mmax -0.003 0.11]);
yticks([0 0.025 0.05 0.075 0.1]);
xticks([0 50 100 150 200 250 300]);
latexylabel('$${\rm Prob}(M=m)$$');
rotateY;
latexxlabel('$$m$$');
%plot([62, 118],[0 0],'b','linewidth',4);   % is changed when generating eps-figures
fill([62,62,118,118,62],[0 -0.002 -0.002 0 0],'b');
hold off
%
figure(2)
axis([0 mmax 0 1.1]);
yticks([0 0.2 0.4 0.6 0.8 1]);
xticks([0 50 100 150 200 250 300]);
latexylabel('$${\rm Prob}(M \le m)$$');
rotateY;
latexxlabel('$$m$$');
hold off


%%  Distribution of random graphs with 3 vertices
%   Theory: Binomial distribution
figure(7)
subplot(3,3,[1,2,4,5]);
N1=3;
mmax=(N1-1)*N1/2;
p=0.4;
%  draw the Binomial distribution
Prob=[];
Sum=[];
Prob(1)=(1-p)^mmax;
for k1=1:mmax
    Prob(k1+1)=nchoosek(mmax,k1)*p^k1*(1-p)^(mmax-k1);     
                     % nchoosek determines the binomial coefficient
end
for k1=0:mmax
    fill([k1-0.25 k1-0.25 k1+0.25 k1+0.25 k1-0.25],[0 Prob(k1+1) Prob(k1+1) 0 0],[0.73 0.91 1]);
    hold on
end
plot(0:mmax, Prob, 'bo',...
            'MarkerEdgeColor','b',...
            'MarkerFaceColor','b');
% %   expected value for the number of edges
% Em=p*mmax;
% plot([Em Em],[0 max(Prob)],'b--');
axis([-0.5 3.5 -0.01 0.5]);
yticks([0 0.1 0.2 0.3 0.4 0.5]);
xticks([0 1 2 3]);
latexylabel('$${\rm Prob}(M=m)$$');
rotateY;
latexxlabel('$$m$$');
hold off

%%  Characteristic path length of random graphs
%   in dependence upon the probability p
%   from  NSC: RandomGraphDemo.m
%
figure(3)
subplot(3,1,[1,2]);
N1=40;  % 40
kmax=700;
Prob=[0.07:0.015:0.3 0.35:0.05:1];
Pathlength=[];            % characteristic path length
                    % in dependence upon the probability p
kp=0;
%RandomGraphs2=rng;                         % save s to restore the generator setting
%save('RandomGraphs2', 'RandomGraphs2');
load('RandomGraphs2');       % Seed (initialisation) of the number generator
rng(RandomGraphs2);
%
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 RandomGraphs:}$$ 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.63,'$$\tilde{p}=0.22$$');
hold off

%%  Characteristic path length of random graphs
%   in dependence upon the vertex number N1
%   from  NSC: RandomGraphDemo.m
%
figure(4)
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 RandomGraphs:}$$ Diameter of random graphs');
latextext(800, 3.5, ['$$\bar{d}=\;$$', num2str(dbar)]);
hold off


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%   Connectivity of random graphs
%
%  determine the lower bound of the connection probability p
%  such that the undirected random graph is probable to be connected
%  demonstrates a phase transition
%  from RandomGraphDemo.m  for NCS
%
figure(5)
subplot(4,1,[1,2])
N1=20;
kmax=450;
Prob=[0.02:0.005:0.39 0.4:0.1:1];
ConnectSum=[];    % average number of connected vertex sets 
                  % in dependence upon the probability p
%RandomGraphs1=rng;                         % save s to restore the generator setting
%save('RandomGraphs1', 'RandomGraphs1');
load('RandomGraphs1');       % Seed (initialisation) of the number generator
rng(RandomGraphs1);
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 RandomGraphs:}$$ Connectivity of random graphs');
hold on
%
%  Determine the probability of the graph to be connected
%
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
subplot(4,1,[1,2])
N1=45;
kmax=150;
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
%
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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Effective resistance of a random electrical network
%
%  from RandUQNet.m
% 
figure(6)
subplot(3,1,[1,2]);
N1=100;                  % Anzahl der Knoten
ReffSum=[];
kmax=300;               % Anzahl der untersuchten Netzwerke für dieselbe Wahrscheinlichkeit
kN=0;
%RandomGraphs3=rng;                         % save s to restore the generator setting
%save('RandomGraphs3', 'RandomGraphs3');
load('RandomGraphs3');       % Seed (initialisation) of the number generator
rng(RandomGraphs3);
Prob=[0.005:0.01:0.6, 1];    % Verbindungswahrscheinlichkeit
for p=Prob
    p                %  Kontrolle des Berechnungsfortschritts
    ReffSum=0;
    k2=0;
    for k=1:kmax
        %  from RandomGraphDemo.m:
        %  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
        Lground=L(1:N1-1,1:N1-1);
        L22=Lground(2:N1-1,2:N1-1); 
        L11=Lground(1,1);
        B=Lground(1,2:N1-1);
         L21=Lground(2:N1-1,1);
        if det(L22)~=0
           k2=k2+1;
           ReffSum=ReffSum+1/(L11-B*inv(L22)*L21);
        end
    end
kN=kN+1;
Reff(kN)=ReffSum/k2;
end
plot(Prob, Reff,'b');
latexxlabel('$$p$$');
latexylabel('$$R_{\rm eff}$$');
rotateY;
xticks([0 0.2 0.4 0.6 0.8 1]);
yticks([0 0.1 0.2 0.3 0.4])
axis([0 1 0 0.3]);
latextitle('$${\tt RandomGraphs:}$$ Resistance of a random network');
hold off
%
%final value for p \to 1
Reff=1/(N1-1-ones(1,N1-2)*inv(N1*eye(N1-2)-ones(N1-2,N1-2))*ones(N1-2,1));

%%
epsfigc('RandomGraphs');