AA_BearRand_master_file.m

## AA_BearRand_master_file.m
%   %   %   %   %   %   %   %   %   %   %   %   %   %   %   %   %   %   %
%                 CODE FOR BEAR RAND (2016) PNAS
%
% This MATLAB script runs the low mutation limit (steady state) calculations
% shown in Figures 2a and S1-S3, as well as the agent-based simulations
% for higher mutation rates shown in Figure S3, and plots the results.
%
% This script calls the following functions:
%
% calculate_payoff(strategy1,strategy2,b,c,d,p)
%       - returns the payoff of strategy1 against strategy2 given
%         the specified values of b, c, d, and p
%
% steady_state_calc(payoffTable,N,w,a)
%       - returns the steady state distribution of the population in
%         the limit of low mutation, where payoffTable is an MxM matrix
%         with payoffTable(i,j) being the payoff of strategy i against
%         strategy j; N is the population size; w is the selection
%         strength; and a is the amount of assortment
%
% agent_based_sim(N,b,c,d,p,w,mu,nGens,a)
%       - returns the average value of each strategy parameter for
%         an agent based simulation using the specified parameter
%         values, run with a mutation rate of mu and last nGens
%         generations.
%   %   %   %   %   %   %   %   %   %   %   %   %   %   %   %   %   %   %

%%%%%%%%%%%%%%%%%%%%%%%
%% DEFINE PARAMETERS %%
%%%%%%%%%%%%%%%%%%%%%%%
clear all

% Parameters used from main text, Figure 2
N = 50; b = 4; c = 1; d = 1; w = 6; a = 0;

% Agent sim parameters from Figure S3a in Supplementary Information
mu = .05; nGens = 1e7; % -> 10mil generations

% We loop through values of p (repeated game probability)
ps_steadyState = .01:.01:.99;
% Because the agent based simulations take much longer, we use a coarser
% grain
ps_agentSim = .1:.1:.9;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% STEADY STATE CALCULATION %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Generate Strategy Space %%%

% For the steady state calculation, we must generate a discrete strategy
% set. To do so:
%   sR, s1, and sI are all restricted to being either 0 or 1;
%   we discretize the deliberation threshold (T) in intervals of d/10

% generate all possible combinations of the first three strat parameters (sR,s1,sI)
% set to be 0 or 1
[x1,x2,x3] = ind2sub([2,2,2], find(zeros(2,2,2)==0));
strats = [x1 x2 x3] - 1;

% we now combine the above with the 11 possible deliberation thresholds,
% running from 0 to d in increments of d/10
Tvals = 0:d/10:d;
strats = repmat(strats,length(Tvals),1);
Tvals = sort(repmat(Tvals,1,length(strats)/length(Tvals)));
strats = [strats Tvals'];

% ... to create a total of 2*2*2*11 (88) possible strategies
M = size(strats,1); % define M to be number of strats


%%% Main Steady State Loop %%%

% initialize matrix to store steady state distributions for each p value
ssDists = zeros(M,length(ps_steadyState));

% initialize matrix to store average strat values for each p value
avgStrats = zeros(length(ps_steadyState),4);

for p=ps_steadyState
    p % display p value that's running

    % initialize payoff table
    payoffTable=zeros(M);

    % create payoff table
    for i=1:M
        for j=1:M
            payoffTable(i,j)=calculate_payoff(strats(i,:),strats(j,:),b,c,d,p);
        end
    end

    % run steady state using this payoff table and store result in ssDists
    ssDists(:,p==ps_steadyState)=steady_state_calc(payoffTable,N,w,a);

    % use ssDists to get average strategy values by taking a weighted average of
    % the strat values based on the distribution in ssDists
    avgStrats(p==ps_steadyState,:) = sum(repmat(ssDists(:,p==ps_steadyState),1,4)...
        .*strats);
end


%%% Plot Results %%%
figure
plot(ps_steadyState,avgStrats); axis([0 1 0 1])
legend('sR','s1','sI','T');
xlabel('Probability of Repeated Game p')
ylabel('Average Value')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% AGENT-BASED SIMULATION %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Main Loop %%%

% initialize matrix to store average strat values for each p value
avgStrats = zeros(length(ps_agentSim),4);

for p=ps_agentSim
    p % display p value that's running

    % run agent sim for select p value
    avgStrats(p==ps_agentSim,:) = ...
        agent_based_sim(N,b,c,d,p,w,mu,nGens,a);
end


%%% Plot Results %%%
figure(2)
plot(ps_agentSim,avgStrats,'--o'); axis([0 1 0 1])
legend('sR','s1','sI','T');
xlabel('Probability of Repeated Game p')
ylabel('Average Value')

## agent_based_sim.m
function [avgStrat,params]=agent_based_sim(N,b,c,d,p,w,mu,nGens,a,popInit)
% Agent based evolutionary dynamics simulations using the Moran process with
% exponential fitness
%
% INPUTS
% N = population size
% w = selection strength
% mu = mutation rate
% nGens = number of generations per simulation
% popInit = Nx4 vector of how you want the population initialized (if this
%           is left off, random initial conditions are generated)
% b,c,d,p are the game parameters used in the payoff function
% a = assortment (probability of interacting with a player having the same
%                 strategy as you, rather than a randomly selected player)
%
% OUTPUTS
% avgStrat = 1x4 vector with cumulative average value of the 4 strategy parameters
%           [sR,s1,sI,T]
% params = matrix of parameters used (optional)v


%%%%%%%
% Initialize population - Nx4 matrix, where
    % strat consists of [sR,s1,sI,T],
    % with sR,s1,sI ranging from 0 to 1 and T ranging from 0 to d
    % (see calculate_payoff.m for details)

% If fewer than 9 inputs were passed to the fn (ie popInit was left
% off) start from a random initial population
if nargin<10
    pop=rand(N,4); % pick random values in [0,1] for each param
    pop(:,4)=pop(:,4)*d; % multiply T by d so it's in [0,d]
    params=[N,b,c,d,p,w,mu,nGens,a];

% Otherwise, use the initial condition you were passed
else
    pop=popInit;
    params=[N,b,c,d,p,w,mu,nGens,a,popInit];

    if N~=length(pop(:,1)) % just in case N and the size of popInit are different, let the user know
                           % and then use the N implied by the initial condition
        disp('N and size of initial pop vector dont match!')
        pause
        N=length(pop(:,1));
    end
end

% Initialize the matrix that will store the payoffs that result from
% each pair of agents interacting
pairwisePayoffs=zeros(N);
for i=1:N
    for j=i:N
        pairwisePayoffs(i,j)=calculate_payoff(pop(i,:),pop(j,:),b,c,d,p);
        pairwisePayoffs(j,i)=calculate_payoff(pop(j,:),pop(i,:),b,c,d,p);
    end
end
% (Leaves diagonal as 0s, since agents dont interact with themselves)

% Initialize the average strategy vector for this simulation run
avgStrat=zeros(1,4);

%%%%%%
% Main simulation loop
for gen=1:nGens

    % Display every 100,000th generation in command window
    if mod(gen,100000)==0
        gen
    end

    % To make the code run faster, we only want to update
    % pairwisePayoffs in a given generation if the population make-up
    % changed. So this variable will keep track of whether, in the
    % curent generation, a change occurred. We initialize it as 0 (ie
    % change did not occur) and then may change it to 1 below if a
    % change is made to pop
    changeOccurred=0;

    % Did a mutation happen?
    if rand<mu
    % Mutation occurs:

        % Pick a learner
        learner=randperm(N,1);

        % Pick a random strategy for her
        pop(learner,:)=rand(1,4); % pick random values in [0,1] for each param
        pop(learner,4)=pop(learner,4)*d;  % multiply T by d so it's in [0,d]

        % Indicate that a change occurred to the population
        changeOccurred=1;

    else
    % No mutation:

        % Pick a random learner
        learner=randi(N);

        % Pick someone to imitate, proportional to fitness
        % (which is exp(w*(a * payoff against own strategy + (1-a)*(avg payoff against all other strategies))))
        teacher=pickOneProportionally(...
            exp(w*(...
                  a*diag(pairwisePayoffs)...   % a * (diagonal of pairwisePayoffs)
                 +(1-a)*sum(pairwisePayoffs.*~eye(N),2)/(N-1))));  %(1-a) * (row-sum of all non-diagonal elements of pairwisePayoff, divided by (N-1) to make it an average)

        % Check if a change will happen
        changeOccurred=(pop(learner,:)~=pop(teacher,:));

        % Update learner's strategy
        pop(learner,:)=pop(teacher,:);

    end

    % Now that we are done with potential mutation/learning, we want to
    % update pairwisePayoffs, but only if pop may have changed (to save
    % computation time)
    if changeOccurred==1
        % go thru each agent, and update the payoff that agent gets
        % against learner, as well as payoff learner gets against that
        % agent
       for i=1:N
            pairwisePayoffs(i,learner)=calculate_payoff(pop(i,:),pop(learner,:),b,c,d,p);
            pairwisePayoffs(learner,i)=calculate_payoff(pop(learner,:),pop(i,:),b,c,d,p);
       end
    end

    % add this generation to avgStrat
    avgStrat=avgStrat+mean(pop);

end
% Now we are done with the main simulation loop - this simulation run
% is over.

% convert avgStrat to an average from sum of averages across all gens
avgStrat=avgStrat/nGens;

end


function [index,value]=pickOneProportionally(input)
% 'input' should be an Nx1 vector of values.
% this fn then picks an entry proportional to value, and
% returns that entry's index and value.

% % set NaNs to 0 - NaNs should never occur, but just
% badIdx=find(isnan(input));
% input([badIdx badIdx+length(input)/2])=0;

% first, check if input is all 0s, in which case you pick one at random
if sum(input)==0
    index=randi(length(input));
    value=0;

else
    % otherwise:

    % normalize values so they sum to 1
    inputN=input/sum(input);

    % make a new vector 'scaled' is the CDF of input
    % (i.e. that effectively lines up the values along the number line from 0 to 1)
    % for example, say that you had
    % input=[4 2 1 3 10]
    % then after normalization you'd have input=[.2 .1 .05 .15 .5]
    % and the CDF would be [.2 .3 .35 .5 1]
    scaled=zeros(length(inputN),1);
    scaled(1)=inputN(1);
    for i=2:length(inputN)
        scaled(i)=scaled(i-1)+inputN(i);
    end

    % now generate a random number between 0 and 1, and see which bin of scaled
    % it lands in: that is, find the first entry in scaled which is larger than
    % the random number
    index=find((scaled>rand)==1,1);
    value=input(index);

end


end

## calculate_payoff.m
function payoff = calculate_payoff(s1,s2,b,c,d,p)
%% SETUP %%

% s1 and s2 are agent 1's and agent 2's respective strategy profiles
%   where a strategy profile consists of [pR,p1,pH,T]

% STRATEGY PARAMETERS %
% pR is probability of cooperating when agent deliberates and realizes the
% game is repeated/reciprocal
% p1 is probability of cooperating when agent deliberates and realizes the
% game is 1-shot/social dilemma
% pI is probability of cooperating when agent uses intuition (and therefore
% can't discriminate repeated vs. 1-shot games)
% T is (threshold) max cost in [0,d] agent is willing to pay to deliberate

% ENVIRONMENT/GAME PARAMETERS %
% b is benefit of cooperation
% c is cost of cooperation
% d is max possible cost of deliberation (such that the cost of
%   deliberation varies uniformly in the interval [0,d])
% p is probability of repeated game

% OUTPUTS %
% payoff is s1's payoff


%%%%%%%%%%%%%%%%%%%%%%%
%% CALCULATE PAYOFFS %%
%%%%%%%%%%%%%%%%%%%%%%%

% Rename each agent's threshold to T1 and T2, respectively
T1 = s1(4);
T2 = s2(4);

%%%%
% Calculate repeated game payoffs

% Repeated Game Payoff Matrix
repPayoffMat=[0	0
              0 (b-c)];

%   Rows correspond to agent 1's strategy
%       where row 1 = agent 1 defects, row 2 = agent 1 cooperates
%   Columns correspond to agent 2's strategy
%       where col 1 = agent 2 defects, col 2 = agent 2 cooperates
%   Values correspond to agent 1's payoffs
%       e.g., if agent 1 cooperates and agent 2 defects,
%       agent 1 gets a payoff of 0 because repPayoffMat(2,1) = 0

% We consider repeated-game payoffs in 4 situations:
%   rPayDD: both agents deliberate
%   rPayDI: only agent 1 deliberates
%   rPayID: only agent 2 deliberates
%   rPayII: neither agent deliberates

% E.g., to calculate rPayDI,
%   We know agent 1 deliberates and realizes they're in a repeated game =>
%       their prob. of cooperating is pR, i.e., s1(1)
%   We know agent 1 uses their intuitive response =>
%       their prob. of cooperating is pI, i.e., s2(3)
%   Thus, rPayDI is a weighted average of the payoff agent 1 gets when both
%       cooperate (repPayoffMat(2,2)), when only player 1 cooperates
%       (repPayoffMat(2,1)), and so on
%   where the prob. both cooperate is given by s1(1)*s2(3), prob. only
%   player 1 cooperates is given by s1(1)*(1-s2(3)), and so on

rPayDD=s1(1)*s2(1)*repPayoffMat(2,2) + s1(1)*(1-s2(1))*repPayoffMat(2,1) + ...
    (1-s1(1))*s2(1)*repPayoffMat(1,2) + (1-s1(1))*(1-s2(1))*repPayoffMat(1,1);
rPayDI=s1(1)*s2(3)*repPayoffMat(2,2) + s1(1)*(1-s2(3))*repPayoffMat(2,1) + ...
    (1-s1(1))*s2(3)*repPayoffMat(1,2) + (1-s1(1))*(1-s2(3))*repPayoffMat(1,1);
rPayID=s1(3)*s2(1)*repPayoffMat(2,2) + s1(3)*(1-s2(1))*repPayoffMat(2,1) + ...
    (1-s1(3))*s2(1)*repPayoffMat(1,2) + (1-s1(3))*(1-s2(1))*repPayoffMat(1,1);
rPayII=s1(3)*s2(3)*repPayoffMat(2,2) + s1(3)*(1-s2(3))*repPayoffMat(2,1) + ...
    (1-s1(3))*s2(3)*repPayoffMat(1,2) + (1-s1(3))*(1-s2(3))*repPayoffMat(1,1);


%%%
% Calculate 1-shot game payoffs
%   (these work the same way as the repeated game payoffs)

% 1-shot Game Payoff Matrix
onesPayoffMat=[0 b
              -c (b-c)];

oPayDD=s1(2)*s2(2)*onesPayoffMat(2,2) + s1(2)*(1-s2(2))*onesPayoffMat(2,1) + ...
    (1-s1(2))*s2(2)*onesPayoffMat(1,2) + (1-s1(2))*(1-s2(2))*onesPayoffMat(1,1);
oPayDI=s1(2)*s2(3)*onesPayoffMat(2,2) + s1(2)*(1-s2(3))*onesPayoffMat(2,1) + ...
    (1-s1(2))*s2(3)*onesPayoffMat(1,2) + (1-s1(2))*(1-s2(3))*onesPayoffMat(1,1);
oPayID=s1(3)*s2(2)*onesPayoffMat(2,2) + s1(3)*(1-s2(2))*onesPayoffMat(2,1) + ...
    (1-s1(3))*s2(2)*onesPayoffMat(1,2) + (1-s1(3))*(1-s2(2))*onesPayoffMat(1,1);
oPayII=s1(3)*s2(3)*onesPayoffMat(2,2) + s1(3)*(1-s2(3))*onesPayoffMat(2,1) + ...
    (1-s1(3))*s2(3)*onesPayoffMat(1,2) + (1-s1(3))*(1-s2(3))*onesPayoffMat(1,1);


%%%
% Calculate agent 1's final payoff based on
%   Prob. of repeated game (p) vs. 1-shot game (1-p)
%   Prob. agents 1 and 2 deliberate (T1/d and T2/d) vs. use intuition (1-T1/d and 1-T2/d)
%   Average costs agent 1 and 2 pay when deliberating (T1/2 and T2/2)
%       (because the cost of deliberation varies uniformly, an agent who's
%       willing to pay max cost of deliberation T pays an average cost of
%       T/2 to deliberate)

% We take a weighted avg. of agent 1's payoff in 4 situations:
%   1. payoff when both agents deliberate = (p*rPayDD + (1-p)*oPayDD - T1/2)
%       this happens with probability (T1/d)*(T2/d)
%   2. payoff when only agent 1 deliberates = (p*rPayDI + (1-p)*oPayDI - T1/2)
%       this happens with probability (T1/d)*(1 - T2/d)
%   ... and so on

payoff = (T1/d)*(T2/d)*(p*rPayDD + (1-p)*oPayDD - T1/2) + ...
    (T1/d)*(1 - T2/d)*(p*rPayDI + (1-p)*oPayDI - T1/2) + ...
    (1 - T1/d)*(T2/d)*(p*rPayID + (1-p)*oPayID) + ...
    (1 - T1/d)*(1 - T2/d)*(p*rPayII + (1-p)*oPayII);


end

## steady_state_calc.m
function [ssDist,transitionPs]=steady_state_calc(payoffTable,N,w,a)
% Calculate the steady state distribution of the Moran process with
% exponential fitness, assuming the limit of low mutation
%
% INPUTS
% payoffTable is an MxM matrix of payoffs, where M is the number of strategies
% N is the population size
% w is the selection strength (0=neutral drift, ->Inf=deterministic)
% a is assortment (probability of interacting with a player having the same
%                 strategy as you, rather than a randomly selected player)
%
% OUTPUTS
% ssDist is an Mx1 vector with the steady state frequency of each strategy
% transitionPs is an MxM vector with transition probilities between
%              homogeneous states of the population

numStrats=length(payoffTable);
transitionPs=zeros(numStrats);

% To do the steady state calc, we need to assemble a transition matrix
% betewen homogeneous states of the population
for s1=1:numStrats
    for s2=1:numStrats

        if s1~=s2

        % for each combination of strategies s1 and s2, the transition prob
        % going from s1 to s2 is the product of the probability a s1 mutant
        % arises, and the probability that mutant goes to fixation

                            %prob s1 mutant arises           %prob that s1 goes to fixation
                            %(assuming uniform mutation)
            transitionPs(s1,s2)=(1/(numStrats-1))            *calc_invasionprob(payoffTable,N,w,s1,s2,a);

        end
    end
end

% The diagonal entries of the transition matrix are the probability you
% stay in the given state - which is just 1 minus all the transition
% probabilities out of that state
for i=1:numStrats
    transitionPs(i,i)=1-sum(transitionPs(:,i));
end

% Now that you've got the transition matrix, find the eigenvector with
% eigenvalue 1, and that's your steady state distribution!
[vec,val]=eig(transitionPs);
[~, idxofEigWithVal1]=min(abs(diag(val)-1));  % For some reason, sometimes Matlab gets the eigenvalue wrong s/t its not literally 1, just very close

% normalize the eigenvector corresponding to eigenvalue so that it sums to 1, and you're done!
ssDist=vec(:,idxofEigWithVal1)./sum(vec(:,idxofEigWithVal1));

end


function rhoBA=calc_invasionprob(payoffTable,N,w,sA,sB,a)
% here we calculate rhoAB, the probability that a sA mutant goes to fixation in
% a sB population

% this probability is the inverse of (1 + a sum of products) so I
% start by intitalizing rhoInv to 1, and then add on the sum of series of
% products, and then take the inverse
rhoInv=1;

% loop through each possible number of steps the sA mutant would need to go through in order to take over
% (this ranges from 1, if there are already N-1 sA players,
% to N-1, if there is only the 1 initial sA mutant)
for k=1:N-1

    % same thing - initialize the eventual product of various terms as 1,
    % then multiply it by stuff
    thisProd=1;

    % loop through each possible make-up of the population (indexed by k=the
    % number of sA players) that would need to be traversed to get from
    % there being only 1 sA player to there being k sA players
    for i=1:k

        % for each population make-up, multiply by the ratio of the fitness
        % of an sA player in the current pop to the fitness of an sB player
        % in the current pop (where there are i sA players)

        % fitness = exp( w * expectedpayoff)
        % expectedpayoff = a * payoff against your own strategy
        %                  + (1-a) *( (payoff against sB * fraction of other players who are sB)
        %                            +(payoff against sA * fraction of other players who are sA)
        %                           )

        fi=exp(w* (a*payoffTable(sA,sA)+(1-a)*(payoffTable(sA,sB)*(N-i)/(N-1)   + payoffTable(sA,sA)*(i-1)/(N-1))));
        gi=exp(w* (a*payoffTable(sB,sB)+(1-a)*(payoffTable(sB,sB)*(N-i-1)/(N-1) + payoffTable(sB,sA)*(i)/(N-1)  )));

        % multiple the current state of the product by the ratio for this value of i
        thisProd=thisProd*(gi/fi);
    end
    % now that we've finished with the product for k, add it onto the sum
    rhoInv=rhoInv+thisProd;
end

% now that we've finished with the sum, take the inverse
rhoBA=1/rhoInv;

end