% Reinforcement learning for a costly signaling game in which there are
% two types. Code written by Carl Veller, 2013. 

% This code calculates timetrends for a population of size pop_size over 
% number of generations given by num_gens. Mutation rate is given by u. 
% The cost of signal to hi type and lo type, respectively, are given 
% by c_lo and c_hi. 


clc
clear all

tic

% Determine number of generations and iterations.
num_gens = 2000;
num_iters = 1;

% Initialize hi and lo cost vectors. 
hi_type_costs = zeros(1, num_iters);
lo_type_costs = zeros(1, num_iters);

% Determine frequency of low type. 
q = 2/3;

% Determine population size.
pop_size = 100;


num_receivers = round(pop_size/2);
num_hi = round(pop_size*(1-q)/2);
num_lo = round(pop_size*(q)/2);

% Determine selection strength.
alpha = .5;
gamma = 2;

% Determine mutation rate. 
u = 0;


% Payoff to sender of signal, for being received. 
payofftosender = 5;

% Costs of signals by type.
c_lo = [0, 3, 6, 9];
c_hi = [0, 1, 2, 3];

[r signals] = size(c_hi);

% Determine payoffs to receivers. 
% (Left corresponds to lo, right corresponds to hi.)
firm_pi = [0, 10];

% Possible action profiles for sender. 
action = [0 0 0 0;
    0 0 0 1;
    0 0 1 0;
    0 1 0 0;
    1 0 0 0;
    0 0 1 1;
    0 1 0 1;
    1 0 0 1;
    0 1 1 0;
    1 0 1 0;
    1 1 0 0;
    0 1 1 1;
    1 0 1 1;
    1 1 0 1;
    1 1 1 0;
    1 1 1 1]; 


for iteration=1:1:num_iters
    time=0;
    
    % Initialize arrays...

    % ...for each strategy
    value_hi = zeros(num_hi, signals, num_gens);
    value_lo = zeros(num_lo, signals, num_gens);
    value_receivers = zeros(num_receivers, 2^signals, num_gens); 
    
    % ...for probabilities of playing each strategy before possible mutation. 
    prob_hi = zeros(num_hi, signals, num_gens);
    prob_lo = zeros(num_lo, signals, num_gens);
    prob_receiver = zeros(num_receivers, 2^signals, num_gens); 

    % ...for number of strategies.
    strats_hi = zeros(num_hi, signals, num_gens);
    strats_lo = zeros(num_lo, signals, num_gens); 
    strats_receivers = zeros(num_receivers, 2^signals, num_gens); 
    actions_receivers = zeros(num_receivers, signals, num_gens);
    
    % ...for 
    payoff_hi = zeros(num_hi, 1, num_gens);
    payoff_lo = zeros(num_lo, 1, num_gens);
    payoff_receiver = zeros(num_receivers, 1, num_gens);
    

    % Determine random initial values.
    value_hi(:,:,1) = rand(num_hi, signals, 1);
    value_lo(:,:,1) = rand(num_lo, signals, 1);
    value_receivers(:,:,1) = rand(num_receivers, 2^signals, 1); 

   
    
    for t = 1:1:num_gens
        
       % Generate strategies from value functions.
                
        exp_value_hi = exp(gamma*value_hi(:,:,t));
        prob_hi(:,:,t) = (exp_value_hi'./(ones(size(exp_value_hi',1),1)*sum(exp_value_hi')))';
        
        for k = 1:num_hi
            choose_strat = randsample(signals,1,true,prob_hi(k,:,t));
            % Do we mutate?
            mutate1 = rand; 
            % If mutate, mutate to this.
            mutate2 = floor(signals*rand) + 1; 
            final_strat = choose_strat*(mutate1>u) + mutate2*(mutate1<u);
            final_strat = choose_strat;
            strats_hi(k,final_strat,t) = 1;
        end
        
        exp_value_lo = exp(gamma*value_lo(:,:,t));
        prob_lo(:,:,t) = (exp_value_lo'./(ones(size(exp_value_lo',1),1)*sum(exp_value_lo')))';
        
        for k = 1:num_lo
            choose_strat = randsample(signals,1,true,prob_lo(k,:,t));
            % Do we mutate? 
            mutate1 = rand; 
            % If mutate, mutate to this.
            mutate2 = floor(signals*rand) + 1; 
            final_strat = choose_strat*(mutate1>u) + mutate2*(mutate1<u);
            final_strat = choose_strat;
            strats_lo(k,final_strat,t) = 1;
        end
        
        exp_value_receivers = exp(gamma*value_receivers(:,:,t));
        prob_receiver(:,:,t) = (exp_value_receivers'./(ones(size(exp_value_receivers',1),1)*sum(exp_value_receivers')))';
        
        for k = 1:num_receivers
            choose_strat = randsample(2^signals,1,true,prob_receiver(k,:,t));
            % Do we mutate?
            mutate1 = rand; 
            % If mutate, mutate to this. 
            mutate2 = floor((2^signals)*rand) + 1; 
            final_strat = choose_strat*(mutate1>u) + mutate2*(mutate1<u);
            final_strat = choose_strat;
            strats_receivers(k,final_strat,t) = 1;            
            actions_receivers(k,:,t) = action(final_strat,:);
        end

        % Calculate payoffs for lo types and hi types.
        choosing_hi = strats_hi(:,:,t)*(actions_receivers(:,:,t))';
        payoff_hi(:,:,t) = payofftosender*(sum(choosing_hi')')/num_receivers - strats_hi(:,:,t)*c_hi';
        
        choosing_lo = strats_lo(:,:,t)*(actions_receivers(:,:,t))';
        payoff_lo(:,:,t) = payofftosender*(sum(choosing_lo')')/num_receivers - strats_lo(:,:,t)*c_lo';
        
        % Determine number of hi types chosen receivers.
        choice_hi1 = actions_receivers(:,:,t)*strats_hi(:,:,t)';
        choice_hi = sum(choice_hi1')';
        
        % Determine number of lo types chosen receivers.
        choice_lo1 = actions_receivers(:,:,t)*strats_lo(:,:,t)';
        choice_lo = sum(choice_lo1')'; 
        
        payoff_receiver(:,:,t) = firm_pi(1)*choice_lo + firm_pi(2)*choice_hi - payofftosender*(choice_lo + choice_hi);
        
        
        % Update value functions for hi type.
        for i = 1:num_hi
               
            if t == 1
                value_hi(i,:,t+1) = value_hi(i,:,t);
            else            
                value_hi(i,:,t+1) = value_hi(i,:,t) + alpha*(payoff_hi(i,:,t) - value_hi(i,:,t)*strats_hi(i,:,t)')*strats_hi(i,:,t);
            end
             
        end
        
        % Update value function for lo type. 
        for i = 1:num_lo
               
            if t == 1
                value_lo(i,:,t+1) = value_lo(i,:,t);
            else            
                value_lo(i,:,t+1) = value_lo(i,:,t) + ...
                    alpha*(payoff_lo(i,:,t) - ...
                    value_lo(i,:,t)*strats_lo(i,:,t)')*strats_lo(i,:,t);
            end
             
        end
        
        % Update value function for receivers. 
        for i = 1:num_receivers
                        
            if t == 1
                value_receivers(i,:,t+1) = value_receivers(i,:,t);
            else            
                value_receivers(i,:,t+1) = value_receivers(i,:,t) +...
                    alpha*(payoff_receiver(i,:,t) - ...
                    value_receivers(i,:,t)*strats_receivers(i,:,t)')*strats_receivers(i,:,t);
            end
             
        end
        
        % Determine total for hi type. 
        total_hi = zeros(1,signals);
        for j = 1:num_hi
          total_hi = total_hi + strats_hi(j,:,t);
        end
        avg_hi(t,:) = total_hi/num_hi;
        avg_cost_hi(t) = avg_hi(t,:)*c_hi';
        
        % Determine total for lo type.
        total_lo = zeros(1,signals);
        for j = 1:num_lo
              total_lo = total_lo + strats_lo(j,:,t);
        end
        avg_lo(t,:) = total_lo/num_lo;
        avg_cost_lo(t) = avg_lo(t,:)*c_lo';
        
        % Determine total for receivers.
        total_receivers = zeros(1,signals);
        for j = 1:num_receivers
               total_receivers = total_receivers + actions_receivers(j,:,t);
        end
        avg_receivers(t,:) = total_receivers/num_receivers;
        
    end
    
    end_avg_cost(iteration) = avg_cost_hi(num_gens);
    total_avg_cost(iteration) = mean(avg_cost_hi);
    
end

% Plot average cost to high type and average cost to low type
plot(1:t,avg_cost_hi, 1:t, avg_cost_lo)
xlabel('Generation')
ylabel('Average cost of signal')

toc





        
        
                                                                    
                                                                    
                                                         
        
        
        
        
        
        
        
        
        

