
%% Comparison of Shift Register and HB in convergence Rate
clear all;
n_iter=300;

N=100;

% % %DummBell
% % undirected

m = floor(N/3);
B1 = full(spdiags([ones(m,1),ones(m,1), ones(m,1)], [0 1 -1], m, m)); B1(1,m)=1; B1(m,1)=1;
B2 =full(spdiags([ones(m,1),ones(m,1), ones(m,1)], [0 1 -1], m, m));B2(1,m)=1; B2(m,1)=1;
B3 =full(spdiags([ones(m,1),ones(m,1), ones(m,1)], [0 1 -1], m, m)); 
A_undirected = zeros(3*m, 3*m); A_undirected(1:m,1:m)= B1; A_undirected(m+1:2*m,m+1:2*m ) = B3; A_undirected(2*m+1:3*m,2*m+1:3*m)=B2; 
A_undirected(m+1, m) = 1; A_undirected(m, m+1) =1;  A_undirected(2*m+1, 2*m) = 1; A_undirected(2*m, 2*m+1)=1;
 

[L,N]=size(A_undirected);
degvec = zeros(N,1);

for i=1:N
    degvec(i, 1)=nnz(A_undirected(i,:)); 
end
degreeMat = diag(degvec);
Q = degreeMat-A_undirected;


z=(1:N)';
%Hb convergence rate
eig(Q); sort(ans); sqrt(ans(end)/ans(2))



%% HB

W = Metropolis(A_undirected,'zero');
eQ=sort(eig(W)); Lmax=eQ(end); lmin=eQ(2);

alpha=4/(sqrt(Lmax)+sqrt(lmin))^2; beta=((sqrt(Lmax)-sqrt(lmin))/(sqrt(Lmax)+sqrt(lmin)))^2;


z=(1:N)';

x_hb=[eye(N,N);eye(N,N)]*z;
xvec_hb = [z];
for k=1:n_iter,
    x_hb=[ (1+beta)*eye(N,N)-alpha*W, -beta*eye(N,N) ; eye(N,N) zeros(N,N) ]*x_hb;
    xvec_hb=[xvec_hb,x_hb(1:N,:)];
    end;

% plot(|x-x^*|^2)
semilogy(sum(((mean(z)*ones(N,1))*ones(1,n_iter+1)-xvec_hb).^2), '-k', 'LineWidth', 2);
hold on;
%%  SHift register


W = Metropolis(A_undirected, 'one');


mu = max(eig(W-ones(N,1)*ones(N,1)'));
zeta = 2/(1+sqrt(1-mu^2));


z=(1:N)';
x=[eye(N,N);eye(N,N)]*z;
xvec=[z];
yvec=[eye(N,N) zeros(N,N)]*x;
for k=1:n_iter,
    x=[zeta*W (1-zeta)*eye(N,N); eye(N,N) zeros(N,N) ]*x;
    y = [eye(N,N),zeros(N,N)]*x;
    xvec=[xvec,x(1:N,:)];
    yvec = [yvec,y];
end;
semilogy(sum(((mean(z)*ones(N,1))*ones(1,n_iter+1)-xvec).^2),  '-.k', 'LineWidth', 2);


hold on;
%% Nesterov

w_metro = Metropolis(A_undirected,'zero');

w = w_metro;
eigs=sort(eig(w)); 
Lmax=max(eigs); lmin=eigs(2);
alpha=1/Lmax; beta=(sqrt(Lmax)-sqrt(lmin))/(sqrt(Lmax)+sqrt(lmin));
a = alpha;
b = beta;
x0 = (1:N)';
x = x0;
x_old = x;
x_vec = [x];
for i = 1:n_iter,
    x_next = (1+b)*(eye(N,N)-a*w)*x - b*(eye(N,N)-a*w)*x_old;
    x_vec = [x_vec, x_next];
    x_old = x;
    x = x_next;
end;
% % leg = [leg;'HB+Metro'];
semilogy(sum((x_vec - ((mean(x0)*ones(N,1))*ones(1,n_iter+1))).^2),'.k', 'LineWidth', 2);
hold on;
%% Standard

W_mh=Metropolis(A_undirected,'one');
z=(1:N)';
x_mh=z;
xvec_mh=[x_mh];
for k=1:n_iter,
    x_mh=W_mh*x_mh;
    xvec_mh=[xvec_mh,x_mh];
end;
semilogy(sum((xvec_mh-(mean(z)*ones(N,1))*ones(1,n_iter+1)).^2), '--k', 'LineWidth',2);
hold on;

%%
legend('multi-step', 'shift-register','Nesterov','standard consensus',1);
% 
% xlabel('t');
%  ylabel('f(x(t))-f^{*}');
hx=xlabel('t');
hy=ylabel('|| x(t)-x^* ||^2');
font_name='times';
% font_name='Candara';

set(hx, 'FontName', font_name, 'FontWeight', 'Bold', 'FontSize', 14);
set(hy, 'FontName', font_name, 'FontWeight', 'Bold', 'FontSize', 14);
set(gca, 'FontName', font_name, 'FontSize', 14, 'FontWeight', 'Demi', 'LineWidth', 2);

set(gca, 'XLim', [0 n_iter]);
