%######################################################################
%%% Run of average consensus algorithm with perfect communication                   
%######################################################################
clear all;
s=[1 1 2 2 3 3 3 4 4 4 5 6 7 8 8 9 9 10  ];%source vertices
t=[2 3 4 5 7 10 5 5 6 7 6 7 5 1 2 10 7 8];%destination vertices
w=[1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]; %weights of edges
G = graph(s,t,w); %create graph G with parameters s, t, w
iterations = 30;

for i= 1:length(s) %A is adjacency matrix; A(i,j)=1 <=> exists edge (i,j)
   A(s(i),t(i))=1;
   A(t(i),s(i))=1; 
end

checkA=isequal(A, adjacency(G)); %chack adjacency using matlab function

M=0; %M is the highest degree of vertex
for i=1:size(A)
    sumRadek=sum(A(i,:))
    if sumRadek>=M
        M=sumRadek;        
    end
end
for i = 1:size(A)
    D(i,i)=sum(A(i,:))%D is degree matrix
end
L=D-A; %Laplacian matrix
checkLaplacian=isequal(laplacian(G),L); %check Laplacian using matlab function
for i= 1:length(s) %Incidence matrix
    E(s(i),i)=1;
    E(t(i),i)=-1;
end
L2=E*transpose(E); %another way to compute Laplacian

checkLaplacian2=isequal(laplacian(G),L2); %check Laplacian using matlab function
I=L*ones(max(s)); %check that Laplacian rows sum up to 1
D = eig(L,'matrix');%diagonal matrix of eigenvalues of laplacian
alpha =1/(2*M); %choose parameter alpha
L_eig=eig(L);
alpha=2/(L_eig(2)+L_eig(max(s)))

for i = 1:max(s)
node_initial_value(i)=i; %initialization of values to average
%running_value(i+1,j)=node_value(j);
running_value(1,i)=node_initial_value(i);
running_value_Error(1,i)=-node_initial_value(i)+mean(node_initial_value);
end
for i= 1:max(s)
final_value(i)=mean(node_initial_value); %expected value in node = average
end
node_value=node_initial_value; %inicialization of nodes values
running_value_Error(1,:)=final_value-node_initial_value;
for i=1:iterations+1
discrete_time(i)= i-1;
end
%Perron_ iteration matrix
for i = 1:iterations;
    IT=eye([max(s) max(s)])-alpha*L;
    rand_vektor = 0.02*(wgn(10,1,0))';
       node_value=node_value+ alpha* rand_vektor;
    node_value= node_value* IT;
    for j=1:max(s)
            running_value(i+1,j)=node_value(j);
            running_value_Error(i+1,j)=final_value(j)-node_value(j);
    %node_value(5)=3; uncomment for convergence to v_3 initial value
    end
end
% plot of the results commands ommited

figure;
plot(discrete_time, running_value(:,:));
ttl1=title( 'Run of algorithm  '  )
ttl1=set(ttl1,'Interpreter','latex','FontSize', 15);
xlbl1=xlabel('Iterations [-]');
xlbl1=set(xlbl1,'Interpreter','latex');
ylbl1=ylabel('Value [-]');
ylbl1=set(ylbl1,'Interpreter','latex');
grid;
figure;
grid on
plot(discrete_time, running_value_Error(:,:));
title('shrinking of error');
grid on
figure;
plot(G); %plot graph G





