clear all
hold off
N=7;  % Number of nodes
k=4;  % Number of edges per node in initial regular graph

if rem(N,2)~=0
    Angle=(2*pi)/(N+1);
    for v=0:N
        y(v+1)=sin(v*Angle);
        x(v+1)=cos(v*Angle);
    end
    %plot(x([1:N+2]),y([1:N+2]))
    %hold on
    plot(x([1:N+1]),y([1:N+1]),'o')
    hold on
    
    % Calculate Regular Graph Measures
    M=zeros(N,N);
    
    for m=1:floor(k/2)
        M=M+diag(ones(N-m,1),m)+diag(ones(m,1),(N-m));
    end
    
    for s=1:N-1
        for t=s+1:N
            if M(s,t)==1
                plot(x([s,t]),y([s,t]))
                hold on
            end
        end
    end
else
    
    Angle=(2*pi)/N;
    for v=0:N
        y(v+1)=sin(v*Angle);
        x(v+1)=cos(v*Angle);
    end
    %plot(x([1:N+2]),y([1:N+2]))
    %hold on
    plot(x([1:N+1]),y([1:N+1]),'o')
    hold on
    
    % Calculate Regular Graph Measures
    M=zeros(N,N);
    
    for m=1:k/2
        M=M+diag(ones(N-m,1),m)+diag(ones(m,1),(N-m));
    end
    
    for s=1:N-1
        for t=s+1:N
            if M(s,t)==1
                plot(x([s,t]),y([s,t]))
                hold on
            end
        end
    end
end
    
    
    
d=zeros(N,N);

for n=1:N
    for r=1:N
        d(n,r)=graphshortestpath(sparse(M),n,r);
    end
end
d(d==Inf)=NaN;
for q=1:N
    geodesics(q)=max(d(q,:));
end

LRegular=mean(geodesics);


% Calculate Random Graph Measures
for p=1:45
    G=M;
    for u=1:20
        for s=1:N-1
            for t=s+1:N
                if G(s,t)>0
                    a=unifSample([0 1],[(1-(0.02*p)) (0.02*p)]);
                    G(s,t)=G(s,t)-a;
                    if a==1
                        c=ceil((N-s)*rand);
                        G(s,s+c)=G(s,s+c)+1;
                    end
                end
            end
        end
        d=zeros(N,N);
        for n=1:N
            for r=1:N
                d(n,r)=graphshortestpath(sparse(G),n,r);
            end
        end
        d(d==Inf)=NaN;
        for q=1:N
            geodesicsR(q)=max(d(q,:));
        end
        MeanOfGeodesicsR(u)=mean(geodesicsR);
    end
    
    LRandom(p)=mean(MeanOfGeodesicsR);
    Ratio(p)=LRandom(p)/LRegular;
    
end

figure
plot(0.01.*[1:2:90],Ratio,'r.')
axis([0 1 0 4])

hold on




% Calculate Random Graph Measures
for p=1:45
    G=M;
    for u=1:100
        for s=1:N-1
            for t=s+1:N
                if G(s,t)>0
                    a=unifSample([0 1],[(1-(0.02*p)) (0.02*p)]);
                    G(s,t)=G(s,t)-a;
                    if a==1
                        c=ceil((N-s)*rand);
                        G(s,s+c)=G(s,s+c)+1;
                    end
                end
            end
        end
        d=zeros(N,N);
        for n=1:N
            for r=1:N
                d(n,r)=graphshortestpath(sparse(G),n,r);
            end
        end
        d(d==Inf)=NaN;
        for q=1:N
            geodesicsR(q)=max(d(q,:));
        end
        MeanOfGeodesicsR(u)=mean(geodesicsR);
    end
    
    LRandom(p)=mean(MeanOfGeodesicsR);
    Ratio(p)=LRandom(p)/LRegular;
    
end

plot(0.01.*[1:2:90],Ratio,'r')
axis([0 1 0 4])
hold on


%view(biograph(M,[1:N],'ShowArrows','off','EdgeType','straight','LayoutType','equilibrium'))
%view(biograph(G,[1:N],'ShowArrows','off','EdgeType','straight','LayoutType','equilibrium'))