% new low lapsly with equality constraint

clear all;

iter_no =1e3;

N=20; l_No=10;
%source rates


% % % % UTILITY
util_function = 'quadratic';
% util_function = 'logarithm';


% % link capacity
c =1e5;

 as = 1e2*ones(N,1);


%link capacity
xs0 = zeros(N,1);
xs = xs0;

%lower -upper bound of source rates;
m = zeros(N,1);
M = c*ones(N,1);

seed = 0.7;

% % % generate a random routing matrix
R = random_route_saturated(l_No,N,seed);


%maximum hop count (distance of the same source)
[l_No, N] = size(R);
Hop_max = 0;
for i=1:N,
    Hop_max = max (Hop_max,nnz(R(:,i)));
end;

% %max Source count (using the same link)
S_max = 0;
for i=1:l_No,
    S_max = max (S_max,nnz(R(i,:)));
end;



if strcmp(util_function , 'quadratic') == 1
delta_2=1;
delta_1=1;

else
delta_1=as(1,1)/((1+c)^2);
delta_2 = 1/delta_1;
end;





L = delta_2*S_max*Hop_max;
l=delta_1*1;




%% Low Lapsley 


xs0 = zeros(N,1);

cl = c*ones(l_No,1);
pl0 = zeros(l_No,1);

gamma = 2/(L+l);




xs0_vec= [xs0];
pl0_vec = [pl0];
for i=1:iter_no,


if strcmp(util_function , 'quadratic') == 1
    xs0 =  max(min(M-R'*pl0,M),m);
else
    xs0= max(min(as./(pl0'*R)' - ones(N,1),M),m);
end;

pl0 = pl0 + gamma*(R*xs0-cl);
pl0_vec=[pl0_vec,pl0];


xs0_vec=[xs0_vec, xs0];
end;

if strcmp(util_function , 'quadratic') == 1
    util_grad = -1/2*(xs-M)'*(xs-M)
else
util_grad = sum(as.*log(ones(N,1)+xs)) 
end;





%% multi-step

alpha = 4/(sqrt(L)+sqrt(l))^2; beta = ((sqrt(L)-sqrt(l))/(sqrt(L)+sqrt(l)))^2;
DM = as./(ones(N,1)+M);
Dm = as./(ones(N,1)+m);

xs1 = xs0;
pl1 = zeros(l_No,1);
pl1_old = pl1;
xs1_vec= [xs1];
pl1_vec = [pl1];


for i=1:iter_no,


%     
if strcmp(util_function , 'quadratic') == 1
    xs1 =  max(min(M-R'*pl1,M),m);
else
    xs1= max(min(as./(pl1'*R)' - ones(N,1),M),m);
end;


pl1_next = pl1 + alpha*(R*xs1-cl)+beta*(pl1-pl1_old);
pl1_old = pl1;
pl1=pl1_next;
pl1_vec=[pl1_vec,pl1];


xs1_vec=[xs1_vec, xs1];


%update alpha beta

end;
if strcmp(util_function , 'quadratic') == 1
    util_Hb = -1/2*(xs1-M)'*(xs1-M);
else
util_Hb = sum(as.*log(ones(N,1)+xs1)); 
end;



%% plots 

semilogy(sum((xs0*ones(1,iter_no+1)-xs0_vec ).^2).^(1/2),'-k', 'LineWidth', 3); hold on;
semilogy(sum((xs1*ones(1,iter_no+1)-xs1_vec ).^2).^(1/2),'--k','LineWidth', 3); hold on;

ylim([1e-10 10*c]);
xlim([0 iter_no]);


legend('low-lapsley', 'multi-step',1);
hx=xlabel('k');
hy=ylabel('||x_s(\mu(k))-x_s^*||_2^2');
font_name='times';

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);

