% spd.m
% spatial Prisoner's Dilemma
%                                                           4/3/3           SG

clear all
close all
warning off
tic

% PARAMETERS

a_11 = 2.6;        % payoff for cooperator against cooperator; assume a_10= 0
a_01 = 2;        % payoff for defector against cooperator;      assume a_00=1
mu= 0.01;      % mutation probability
n = 100;         % system size
T = 500;        % duration

P = ones(n);                        % payoffs
V = [0 0; -1 0; 1 0; 0 1; 0 -1];   % self, above, below,right,left
XX = ones(n);

% INITIAL CONDITIONS

%X = floor(2*rand(n));       % random initial distribution: X=0 (defect) and X=1 (cooperate)
X = zeros(n);                  % no cooperators initially

if 0
	r = 2;                          % radius of initial pocket
	X = zeros(n);
	X(n/2:n/2+r,n/2:n/2+r) = 1;
end

image(X,'EraseMode','background')
set(gca,'DrawMode','fast','XTick',[],'YTick',[])
title('t = 0')
axis square
%M(1) = getframe;


for t=1:T
	
	X1=X(:,[n 1:n-1]);      % states of neighbors left
	X2=X(:,[2:n 1]);         % states of neighbors right
	X3=X([n 1:n-1],:);      % states of neighbors above
	X4=X([2:n 1],:);         % states of neighbors below
	
	B = X1+X2+X3+X4;        % number of cooperators in neighborhood
	
	P = (X==0).*(B*a_01+4-B)+(X==1).*B*a_11;     % payoff 
	
	PP(1,:,:)=P([n 1:n-1],:);      % pay off of neighbors above
	PP(2,:,:)=P([2:n 1],:);         % payoff of neighbors below
	PP(3,:,:)=P(:,[2:n 1]);         % payoff of neighbors right
	PP(4,:,:)=P(:,[n 1:n-1]);      % payoff of neighbors left
	
	for i=1:n
        for j=1:n
            k=randperm(4);               % random order of neighbors
            z=squeeze(PP(:,i,j));         % vector of neighbors payoffs
            p=P(i,j);                           % own payoff
            ind = 1;
            for K=1:4           % if the K(k)th neighbor payoff is larger,accept its state
                if z(k(K))>p  p=z(k(K)); ind=k(K)+1; end
            end
            Next(i,j) = X( mod(i-1+V(ind,1),n)+1, mod(j-1+V(ind,2),n)+1   );
        end
	end
      
	X = Next;
	
	if mu
        XM = X;
        Q = rand(n) < mu;         % sites where mutation occurs
        XM(Q&(X==0)) = 1;          % defector mutates to cooperator
        XM(Q&(X==1)) = 0;          % cooperator mutates to defector
        X = XM;
	end
	
	coop = sum(sum(X))/n^2;
	
	image(XX+X,'EraseMode','background')
	set(gca,'DrawMode','fast','XTick',[],'YTick',[])
	colormap([1 0 1; 0 1 0]);
	title(['t = ', num2str(t),', freq. of coop= ',num2str(coop,2)])
	axis square
	pause(.01)
	
	%M(t+1) = getframe;

end

toc

%movie(M)