clear all; close all; clc;

% Code for problem 2 of Problem set 8. 


% initilization
T = 3;
c = [1;3;2]; % non-constant price needed to get non-trivial answer

%Nvec=[3;10;20;40;60;80;100];
N = 10;


%% use quadratic programming to find the Nash and social welfare optimal
% Nash equilibrium
B1 = ones(N,N) + eye(N,N); % the matrix given in the homework hint
H_p = kron(B1, eye(T,T)); % stacked gradients linear term
f_p = kron(ones(N,1), c); % stacked gradients affine term
% stacked gradients is H_px + f_p, where x is the stacked strategies

LB = zeros(N*T,1); % lower bound is 0
UB = ones(N*T,1); % upper bound is 1
% now equality constraint for each player
Aeq = kron(eye(N,N),ones(1,T));
beq = ones(N,1);
% since the game is potential, we can find the potential function optimizer
[x_NE, P_NE, EXIT_flag]= quadprog(H_p,f_p,[],[],Aeq,beq,LB,UB);




%% using  Social Welfare optimal
B1 = 2*ones(N,N); % the gradients of social welfare is different
H_s = kron(B1, eye(T,T)); % quadratic term in the social welfare
%H_s = 2*ones(T*N, T*N);
f_s = kron(ones(N,1), c);   % affine term in the social welfare
% so the social welfare is x'H_sx + f_s
[x_SW, J_SW]= quadprog(H_s, f_s,[],[],Aeq,beq,LB,UB);

% welfare function evaluatd at NE
J_NE = (x_NE'*H_s*x_NE)/2 + f_s'*x_NE;
J_SW2 = (x_SW'*H_s*x_SW)/2 + f_s'*x_SW


%% using gradient descent algorithm to find Nash
L = (2+N); % Here, Lipschitz constant is the maximum eigenvalue of H_p
mu = 2; % Here, strong monotonicity constant is minimum eigenvalue of H_p
gamma_max = 2*mu/L^2; % max stepsize for convergence (sufficient condition)
K=1000; % number of iterations
ga = 2*gamma_max/3; % pick a stepsize at most gamma_max

x = zeros(T*N,K+1);
x(:,1) = rand(T*N,1); % random initilization
H_proj = eye(N*T, N*T); % the quadratic matrix in projection
for k = 1:K
    y = x(:,k) - ga*(H_p*x(:,k)+f_p); % gradient step
    y_proj = quadprog(H_proj, -y, [],[], Aeq, beq, LB, UB); % projection
    x(:,k+1) = y_proj;

end

% 2-norm distance from the NE
error = norm(x_NE-x(:,end),2) % error in NE estimation

% price of anarchy
PofA = J_NE/J_SW
