%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Competitive and Cooperative Foraging Games: Static One-Stage Case
% By: K. Passino
% Version: 6/17/02
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all

N=2;	% Number of players (can't change)
M=2;	% The number of different types of resources (e.g., nutrients)

% Set the number of cells (locations) in the space (i.e., bins along a line)

Q=21;

% Set the number of different possible values of the decision variables for the players

D1=Q; % Sets that player 1 can make decisions i=1,2,...,D1. Decision corresponds to moving the 
	  % player to position i.
D2=D1;  % Sets that player 2 can make decisions j=1,2,...,D2, meaning placing player 2 at position j.
theta1=1:D1;
theta2=1:D2;

% Set consumption amount. Assume that this same amount is consumed of each type of 
% resource in each cell.

z1=1;	% Sets the amount of resources consumed on one visit to a cell by player 1
z2=1;	% Same, but for player 2

% Define the initial resource map for each type of resource:

for m=1:M
%      r(:,m)=rand(Q,1); % Just a way to initialize
	   r(:,m)=[1 2 3 4 5 6 7 8 9 10 11 10 9 8 7 6 5 4 3 2 1]'; % An example profile (for testing no energy case)
%	   r(:,m)=[1 2 3 4 5 4 3 2 1  2  3  4 5 6 6 6 5 4 3 2 1]'; % An example profile (for testing energy case, two humps, L=1)
	   if m==1
		   r(:,m)=[1 2 3 4 5 4 3 2 1  0  0  0 0 0 0 0 0 0 0 0 0]'; % First resource is close to 0
	   end
	   if m==2
		   r(:,m)=[0 0 0 0 0 0 0 0 1  2  3  4 5 6 6 6 5 4 3 2 1]'; % Second resoure is further away
	   end
end

% Set depletion rate for resource m (independent of cell)

alpha=1*ones(M,1);

% Initialize consumption amounts

C=0*ones(N,M,D1,D2);

% Set resource priorities for each player
p1=ones(M,1);
p2=p1;
p2=[1 2]'; % To make it so that player 2 prefers resource 2

% Next set how much it costs in energy to go get a resource (line thinking of players starting at 0)
% (scales its importance relative to value of resources (so scaling must be correct or it won't move))

we1=0.1;  
we2=0.1;
%we1=0;  % To remove the effect of energy
%we2=0;

% Next, compute how much is consumed of each resource m,  
% if players 1 and 2 make all possible moves (decisions, i and j)

for m=1:M
	for i=1:D1
	for j=1:D2
		if i ~= j
			C(1,m,i,j)=r(i,m)*(1-exp(-alpha(m)*z1)); % Set consumption for forager 1 when forager 2 not there
			C(2,m,i,j)=r(j,m)*(1-exp(-alpha(m)*z2)); % Set consumption for forager 2 when forager 1 not there
		end
		if i==j
			C(1,m,i,j)=0.5*r(i,m)*(1-exp(-alpha(m)*(z1+z2))); % Split return when both present
			C(2,m,i,j)=0.5*r(j,m)*(1-exp(-alpha(m)*(z1+z2))); 
		end
	end
	end
end

% Assume a movement cost that corresponds to each player being at position 0 and having to move to i (j)

% Set weighting matrices on cost of movement (amount of energy to forage):

for i=1:D1
	J1e(i)=we1*i;  % So moving one unit costs we1
end

for j=1:D2
	J2e(j)=we2*j;
end

% Set the cost matrices J1(i,j) and J2(i,j):

for i=1:D1
	for j=1:D2		
			temp1=0;
			temp2=0;
			for m=1:M % Sum up resources consumed for play i,j at location q by both players
				temp1=temp1+p1(m)*C(1,m,i,j);
				temp2=temp2+p2(m)*C(2,m,i,j);
			end
									
		J1(i,j)=-sum(temp1)+J1e(i); % Sums over all cells the consumption for play i,j and the energy 
									  % expended to get it. Note - sign.
		J2(i,j)=-sum(temp2)+J2e(j);
	end			
end			
						
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(1)
clf
bar(1:Q,r,'stacked')
colormap(cool)
xlabel('Cell number')
ylabel('Resource amount')
title('Initial resource distribution')

figure(2)
clf
surf(1:D2,1:D1,J1), colormap(white)
xlabel('Player 2 decision, j')
ylabel('Player 1 decision, i')
zlabel('J_1')
title('Cost function for player 1')

figure(3)
clf
surf(1:D2,1:D1,J2), colormap(white)
xlabel('Player 2 decision, j')
ylabel('Player 1 decision, i')
zlabel('J_2')
title('Cost function for player 2')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute the Nash equilibria:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

flag=0; % Flag for saying if there is no Nash equilibria

for i=1:D1
    for j=1:D2	
		if J1(i,j)<=min(J1(:,j)) & J2(i,j)<=min(J2(i,:)), % Conduct two inequality tests
			display('Nash equilibrium and outcome:') % If satisfied, then diplay solution
			i
			j
			J1(i,j) 
			J2(i,j)
			flag=1; % Indicates that there was one Nash equilibrium (or more)
		end
    end
end

if flag==0
	display('There were no Nash equilibria')
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Minimax strategy
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Compute the minimax strategy (the name "minimax" clearly arises from the computations involved in 
% each player trying to minimize its maximum loss):

% First, compute the security strategy for player 1

[maxvals1,indexmax1]=max(J1'); % This is the max value for each row (note transpose)
[minimaxvalue1,secstratP1]=min(maxvals1) % Display the security value of P1 and its security strategy
                                         % (P1 loses no more than minimaxvalue1, no matter what strategy 
										 % P2 uses)

% Second, compute the security strategy for player 2 (note difference from computation of security
% strategies for matrix games since both payoff matrices are "loss" matrices, and note that computation of the
% minimax strategies is completely independent)

[maxvals2,indexmax2]=max(J2); % This is the min value for each column (note no transpose)
[minimaxvalue2,secstratP2]=min(maxvals2) % Display the security value of P2 and its security strategy

% The outcome of the game will be (with a minimax strategy)

outcome1=J1(secstratP1,secstratP2)
outcome2=J2(secstratP1,secstratP2)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Stackelberg strategy
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Compute the Stackelberg strategy:

% First, find follower reactions (note that this just finds the min loss by P2, so it can
% be that more than one P2 strategy results in this loss):

for i=1:D1
	minJ2(i)=min(J2(i,:)); % Finds the minimum loss of P2 given the leader P1 chooses i
end

P1loss=-inf*ones(1,D1); % Initialization so that it will pick the first min value above it
jstar=0*ones(1,D1);     % Initialization to nonvalid values (valid ones picked next)

for i=1:D1
	for j=1:D2
	if J2(i,j)==minJ2(i)   % Test that j corresponds to a min point for P2, for a given i
		if J1(i,j)>=P1loss(i)  % Tests if it is above any previously stored value for the loss
								% (this finds the security value against all possible P2 rational reactions)
		P1loss(i)=J1(i,j);  % Keeps the value to compare to any other possible P2 reactions later
		jstar(i)=j; % Save the index of the reaction of P2 that results in worst loss for P1 if it uses i
		end
	end
	end
end

% Specify the Stackelberg strategy:

[stackelbergcost,istar]=min(P1loss) % Display the P1 Stackelberg strategy and cost
jstar(istar)                           % Display the P2 follower strategy

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Pareto solutions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% First, start with the Pareto cost defined via scalarization

p=0:0.01:1; % Used for defining the Pareto cost matrices

for k=1:length(p)
	for i=1:D1
		for j=1:D2
			Jp(i,j)=p(k)*J1(i,j)+(1-p(k))*J2(i,j); % Form the Pareto cost matrix
		end
	end
	[temp,indicesmin]=min(Jp); % Find the min elements of each column of Jp, and their indices
	[paretooptimalval(k),paretooptimalstrat2(k)]=min(temp); % Gives min value and its column index
	paretooptimalstrat1(k)=indicesmin(paretooptimalstrat2(k)); % Gives the row index
	J1opt(k)=J1(paretooptimalstrat1(k),paretooptimalstrat2(k));  % Find the cost to J1 for the Pareto-point
	J2opt(k)=J2(paretooptimalstrat1(k),paretooptimalstrat2(k));  % Find the cost to J2 for the Pareto-point
end

figure(4)
clf
subplot(311)
plot(p,paretooptimalval,'ko',p,J1opt,'k-',p,J2opt,'k--')
ylabel('Costs')
title('J_p (o), J_1 (-), and J_2 (--) for Pareto points')
subplot(312)
plot(p,paretooptimalstrat1,'kx')
ylabel('i^*')
title('Decision of player 1 (x)')
subplot(313)
plot(p,paretooptimalstrat2,'k*')
xlabel('Pareto parameter, p')
ylabel('j^*')
title('Decision of player 2 (*)')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute all Pareto solutions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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


%%%%%%%%%%%%%%%%%%% Plot the Pareto points

figure(5)
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('j')
ylabel('i')
title('J_1, J_2, "x" marks a Pareto solution')
hold off


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