%
% The surface plot shows the game solution is the saddle point of the 
% expected payoff function for two strategies game with payoff matrix
% A=[[-2 2]; [4 -3]]
%
clear all
global A W
A=[[-2 2]; [4 -3]];
%A=[[-1 0]; [1 -2]];
%A=[[-2 2]; [0 3]]; %dominent strategy game.
%A=[[-0 0]; [0 3]]; %weakly dominent strategy game.
%A=[[0 -1];[1 -2]];
A=[[2 1];[1 2]]
a=A(1,1);
b=A(1,2);
c=A(2,1);
d=A(2,2);
figure(1)
hold off

[x,y]=meshgrid(0:.01:1);
E=y.*(a*x+c*(1-x))+(1-y).*(b*x+d*(1-x));
meshc(x,y,E);
axis([0 1 0 1 -3 4])
axis([0 1 0 1 1 2])
xlabel('x_1','fontsize',20)
ylabel('y_1','fontsize',20)
zlabel('E','fontsize',20)
hold on

%%
figure(4)
hold off
W=ones(2,2);
[x1,y1]=meshgrid(0:.01:1);

phi11=W(1,1)*max(0*x1,-(1-x1).*((c-a)*y1+(d-b)*(1-y1)));
phi12=W(1,2)*max(0*x1,x1.*((c-a)*y1+(d-b)*(1-y1)));
phi21=W(2,1)*max(0*y1,-(1-y1).*((a-b)*x1+(c-d)*(1-x1)));
phi22=W(2,2)*max(0*y1,y1.*((a-b)*x1+(c-d)*(1-x1)));
V=(phi11.^2+phi12.^2+phi21.^2+phi22.^2)/2;
meshc(x1,y1,V);
axis([0 1 0 1 0 .5])
xlabel('x_1','fontsize',20)
ylabel('y_1','fontsize',20)
zlabel('V','fontsize',20)
hold on

%%
N=5000;
eps=.01;
var0=[0.1867    0.7030];
W=eps*ones(2,2);
%W=eps*rand(2,2);
%var0=[7/11    5/11];
var0=rand(1,2);
var0=[1-0.9891 1-0.0669];
orbit0=var0;
orbit=[];
for i=1:N
    orbit=[orbit;orbit0];
    orbit0=nashmap(orbit0,A,W);
end

%%
offset=.01;
figure(1)
hold on
x=orbit(:,1);
y=orbit(:,2);
E=y.*(a*x+c*(1-x))+(1-y).*(b*x+d*(1-x));
plot3(x,y,E+offset,'k.')
plot3(x(1),y(1),E(1)+offset,'k.','markerfacecolor','w','markersize',30)

%%

figure(4)
hold on
W=ones(2,2);

x1=orbit(:,1);
y1=orbit(:,2);
phi11=W(1,1)*max(0*x1,-(1-x1).*((c-a)*y1+(d-b)*(1-y1)));
phi12=W(1,2)*max(0*x1,x1.*((c-a)*y1+(d-b)*(1-y1)));
phi21=W(2,1)*max(0*y1,-(1-y1).*((a-b)*x1+(c-d)*(1-x1)));
phi22=W(2,2)*max(0*y1,y1.*((a-b)*x1+(c-d)*(1-x1)));
V=(phi11.^2+phi12.^2+phi21.^2+phi22.^2)/2;

plot3(x1,y1,V+offset,'k.')
plot3(x1(1),y1(1),V(1)+offset,'k.','markerfacecolor','w','markersize',30)


%%
figure(2)
hold off
plot(x,y,'bo')
hold on
%plot([7/11],[5/11],'r.','markersize',20)
axis([0 1 0 1])
%return
%%

options=odeset('e','off','BDF','on','AbsTol',1e-8,'RelTol',1e-8);
Tend=50;
numruns=1;
figure(1)
hold on
%W=ones(2,2);
for k=1:numruns
xinit=rand(1,2);
xinit=[0.9891 1-0.0669];
[t,X]=ode15s('nashvectorfield', [0,Tend], xinit, options);

x=X(:,1);
y=X(:,2);
E=y.*(a*x+c*(1-x))+(1-y).*(b*x+d*(1-x));
plot3(x,y,E+offset,'m--','linewidth',4)
plot3(x(1),y(1),E(1)+offset,'m.','markersize',30)
end
view([-25, 40])

oldSettings = fillPage(1, 'margins',[.0 .0 .0 .0]);
%oldSettings = fillPage(1, 'papersize',[8 8]);
%%
figure(4)
hold on
W=ones(2,2);
for k=1:numruns
%xinit=rand(1,2);
%[t,X]=ode15s('nashvectorfield', [0,Tend], xinit, options);

x1=X(:,1);
y1=X(:,2);

phi11=W(1,1)*max(0*x1,-(1-x1).*((c-a)*y1+(d-b)*(1-y1)));
phi12=W(1,2)*max(0*x1,x1.*((c-a)*y1+(d-b)*(1-y1)));
phi21=W(2,1)*max(0*y1,-(1-y1).*((a-b)*x1+(c-d)*(1-x1)));
phi22=W(2,2)*max(0*y1,y1.*((a-b)*x1+(c-d)*(1-x1)));
V=(phi11.^2+phi12.^2+phi21.^2+phi22.^2)/2;

plot3(x1,y1,V+offset,'m--','linewidth',4)
plot3(x1(1),y1(1),V(1)+offset,'m.','markersize',30)
end
view([-25, 40])

% To make high resolution picture, use: oldSettings = fillPage(6000, 'papersize',[10 10]);
% and oldSettings = fillPage(6000, 'margins',[.0 .0 .0 .0]);

oldSettings = fillPage(4, 'margins',[.0 .0 .0 .0]);
%oldSettings = fillPage(4, 'papersize',[8 8]);

return
%%
% Linear Programming method to find NE
%

%A=[[-1 0]; [1 -2]];
%A=[[-2 2]; [4 -3]];
offset=min(min(A));
AA=A-offset;

[r,c]=size(AA);

f=zeros(r+1,1);
f(end)=1;

rone=ones(r,1);
GA=[AA -rone];

b=zeros(r,1);

Aeq=1-f';
beq=1;

lb=0*b;

%options = optimset('LargeScale','off','Simplex','on','Display','iter');
options = optimset('LargeScale','off','Simplex','on');

[optsolu,fval,exitflag,output,lambda] = ...
    linprog(f,GA,b,Aeq,beq,lb,[],[],options);
shadow_price_inequalities=lambda.ineqlin;

plot3([shadow_price_inequalities(1)],[optsolu(1)],[fval+offset],'r.','markersize',30)

%plot3([7/11],[5/11],[35/11-3],'r.','markersize',20)
%%
options=odeset('e','off','BDF','on','AbsTol',1e-16,'RelTol',1e-16);
Tend=500;
xinit=rand(1,2);
xinit=[0.75    0.5001];
[t,X]=ode15s('nashvectorfield', [0,Tend], xinit, options);

x=X(:,1);
y=X(:,2);
figure(2)
hold on
plot(x,y,'k--','linewidth',4)
hold on
plot([shadow_price_inequalities(1)],[optsolu(1)],'r.','markersize',30)


