%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Example of computing all the Pareto solutions for a finite
% game (of course this is using a discrete approximation).
%
% Author: K. Passino
% Version: 4/5/02
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Compute all Pareto solutions for a finite game

clear all


% Define payoff (cost) functions for each player (much coarser)

theta1=-4:0.5:4; % Ok, while we think of it as an infinite game computationally
	              % we of course only study a finite number of points.
m=length(theta1);
theta2=-5:0.5:5;
n=length(theta2);

% A set of cost functions, J1 for player 1 and J2 for player 2
for ii=1:length(theta1)
	for jj=1:length(theta2)
		J1(ii,jj)=-2*(exp( (-(theta1(ii)-2)^2)/5 +(-4*(theta1(ii)*theta2(jj))/20) + ((-(theta2(jj)-3)^2)/2))); 
		J2(ii,jj)=-1*(exp( (-(theta1(ii)-1)^2)/4 + (5*(theta1(ii)*theta2(jj))/10) + ((-(theta2(jj)+1)^2)/2)));
	end
end

PP=ones(size(J1)); % Initialize a matrix that will hold flags indicating if a point is a Pareto Point
					% (initially it indicates that they are all Pareto points)

% Compute the family of Pareto-optimal strategies:

for ii=1:length(theta1)      % These are the loops for the test points theta^*
	for jj=1:length(theta2)

for iii=1:length(theta1)     % These are the loops for the points theta
	for jjj=1:length(theta2)	
		% Perform tests to determine if (ii,jj) is a Pareto point
		if (iii ~=ii & jjj~=jj) &...
			((J1(iii,jjj) <= J1(ii,jj)) & (J2(iii,jjj) <= J2(ii,jj))) &...
			((J1(iii,jjj) < J1(ii,jj)) | (J2(iii,jjj) < J2(ii,jj)))
			PP(ii,jj)=0; % If find one such time that the conditions hold then it is not a Pareto point
		end
	end
end		
		
	end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(1)
clf
colormap(jet)
contour(theta2,theta1,J1,14)
hold on
contour(theta2,theta1,J2,14)
hold on
for ii=1:length(theta1)      % These are the loops for the test points theta^*
	for jj=1:length(theta2)
		if PP(ii,jj)==1
			plot(theta2(jj),theta1(ii),'kx')
			hold on
		end
	end
end
xlabel('\theta^2')
ylabel('\theta^1')
title('J_1, J_2, R_1 (-), R_2 (--), "x" marks a Pareto solution')
hold off


%-------------------------------------
% End of program
%-------------------------------------

