% Individual Learning Model Code - Matlab
% matt.grove@liv.ac.uk
% 2017
% ILFinal.m

% Save the code in Grove_ESM1.txt as an .m function file before running the figure creation sections of this code

clear, clc, close all
tic;

% Setup
ni = 6000; % number of iterations
na = 500; % number of agents
v = 1; % fitness operator variance
mr = 0.05; % mutation sigma
t = (1:ni)'; % time axis
p = 200; % period of wave
st = (1:ni+(p/4))'; % time axis for wave generation

% Parameters
r = [0.1 0.5 0.9]; % reproductive rate
c = [0.1 0.5 0.9]; % cost of IL
s = 0:0.05:1; % learning error (sigma)
a = 0:2.5:50; % amplitude of wave (RoEC)

% Output matrices
[workM,workB] = deal(zeros(ni,3));
[Mphen,Mil,Mfit,Bphen,Bil,Bfit] = deal(zeros(21,21));
[AllM,AllB,AllMFit,AllBFit] = deal(zeros(63,63));
y1 = [43 22 1]; y2 = [63 42 21]; % output indices (y = CoIL)
x1 = [1 22 43]; x2 = [21 42 63]; % output indices (x = RR)

% LOOP 1 – cost of individual learning c(i)
for i = 1:3

% LOOP 2 – reproductive rate r(j)
for j = 1:3

Rn = na*r(j); % reproductive number
S = na-Rn; % number of survivors
bn = ceil(sqrt(S)); % bin number for histograms (for CtM)
 
% LOOP 3 – learning error s(k)
for k = 1:21
 
% LOOP 4 – environmental change a(m)
for m = 1:21

w = a(m)*sawtooth(2*pi*st/p,0.5); % wave generation
w = w(p/4:length(w)-1); % ensures wave starts at zero
 
% Agent setup / reset
CtMagents = zeros(na,3);
CtMagents(1:na,1) = normrnd(0,mr,na,1); % starting phenotype
CtMagents(1:na,2) = 0.1*rand(na,1); % starting proportion individual learning
CtBagents = zeros(na,3);
CtBagents(1:na,1) = normrnd(0,mr,na,1); % starting phenotype
CtBagents(1:na,2) = 0.1*rand(na,1); % starting proportion individual learning
[N,edges] = histcounts(CtMagents(:,1),bn); % CtM setup
[~,I] = max(N);
M = 0.5*(edges(I)+edges(I+1));
B = 0; % CtB setup

for n = 1:ni
CtMagents(:,1) = (1-CtMagents(:,2)).*normrnd(M,s(k),na,1)+CtMagents(:,2).*normrnd(w(n),s(k),na,1); % calculate phenotype
CtBagents(:,1) = (1-CtBagents(:,2)).*normrnd(B,s(k),na,1)+CtBagents(:,2).*normrnd(w(n),s(k),na,1);
CtMagents(:,3) = (1-c(i)*CtMagents(:,2)).*normpdf(CtMagents(:,1),w(n),v); % calculate fitness
CtBagents(:,3) = (1-c(i)*CtBagents(:,2)).*normpdf(CtBagents(:,1),w(n),v);
CtMagents = sortrows(CtMagents,-3); % sort by fitness
CtBagents = sortrows(CtBagents,-3);
%workM(n,1) = median(CtMagents(:,1)); % median phenotype
%workB(n,1) = median(CtBagents(:,1));
workM(n,2) = median(CtMagents(:,2)); % median learning error
workB(n,2) = median(CtBagents(:,2));
workM(n,3) = median(CtMagents(:,3)); % median fitness
workB(n,3) = median(CtBagents(:,3));
Mseed = randsample(S,Rn,true,CtMagents(1:S,3)); % Sample from survivors by FPS
Bseed = randsample(S,Rn,true,CtBagents(1:S,3));
CtMagents(S+1:na,2) = CtMagents(Mseed,2)+normrnd(0,mr,Rn,1); % create & mutate offspring
CtBagents(S+1:na,2) = CtBagents(Bseed,2)+normrnd(0,mr,Rn,1);
Mtb = find(CtMagents(:,2) > 1); % Remove proportion IL > 1
Btb = find(CtBagents(:,2) > 1);
CtMagents(Mtb,2) = 1; % and reset to 1
CtBagents(Btb,2) = 1;
Mtb = find(CtMagents(:,2) < 1); % Remove proportion IL < 0
Btb = find(CtBagents(:,2) < 1);
CtMagents(Mtb,2) = 0; % and reset to 0
CtBagents(Btb,2) = 0;
[N,edges] = histcounts(CtMagents(1:S,1),bn); % find majority phenotype
[~,I] = max(N);
M = 0.5*(edges(I)+edges(I+1));
B = CtBagents(1,1); % find best phenotype
end

% Gather data #1
%Mphen(k,m) = median(workM(1001:ni,1));
Mil(k,m) = median(workM(1001:ni,2));
Mfit(k,m) = median(workM(1001:ni,3));
%Bphen(k,m) = median(workB(1001:ni,1));
Bil(k,m) = median(workB(1001:ni,2));
Bfit(k,m) = median(workB(1001:ni,3));

end
end

% Gather data #2
AllM(y1(i):y2(i),x1(j):x2(j)) = Mil;
AllMFit(y1(i):y2(i),x1(j):x2(j)) = Mfit;
AllB(y1(i):y2(i),x1(j):x2(j)) = Bil;
AllBFit(y1(i):y2(i),x1(j):x2(j)) = Bfit;

end
end

T = toc;

save OP.mat AllM AllMFit AllB AllBFit T

%% Figures
% Proportion of IL; Copy the Best
clear, clc, close all
load OP.mat
a = min(AllB(:)); b = max(AllB(:));
xdata = 0:0.05:1;
zdata1 = AllB(1:21,1:21); zdata2 = AllB(1:21,22:42); zdata3 = AllB(1:21,43:63);
zdata4 = AllB(22:42,1:21); zdata5 = AllB(22:42,22:42); zdata6 = AllB(22:42,43:63);
zdata7 = AllB(43:63,1:21); zdata8 = AllB(43:63,22:42); zdata9 = AllB(43:63,43:63);
Figtype1(a,b,xdata,zdata1,zdata2,zdata3,zdata4,zdata5,zdata6,zdata7,zdata8,zdata9);

%% Proportion of IL; Copy the Majority
a = min(AllM(:)); b = max(AllM(:));
xdata = 0:0.05:1;
zdata1 = AllM(1:21,1:21); zdata2 = AllM(1:21,22:42); zdata3 = AllM(1:21,43:63);
zdata4 = AllM(22:42,1:21); zdata5 = AllM(22:42,22:42); zdata6 = AllM(22:42,43:63);
zdata7 = AllM(43:63,1:21); zdata8 = AllM(43:63,22:42); zdata9 = AllM(43:63,43:63);
Figtype1(a,b,xdata,zdata1,zdata2,zdata3,zdata4,zdata5,zdata6,zdata7,zdata8,zdata9);

%% Proportion of IL; Difference
All = AllB-AllM;
a = min(All(:)); b = max(All(:));
xdata = 0:0.05:1;
zdata1 = All(1:21,1:21); zdata2 = All(1:21,22:42); zdata3 = All(1:21,43:63);
zdata4 = All(22:42,1:21); zdata5 = All(22:42,22:42); zdata6 = All(22:42,43:63);
zdata7 = All(43:63,1:21); zdata8 = All(43:63,22:42); zdata9 = All(43:63,43:63);
Figtype1(a,b,xdata,zdata1,zdata2,zdata3,zdata4,zdata5,zdata6,zdata7,zdata8,zdata9);

%% Fitness; Difference
All = AllBFit - AllMFit;
a = min(All(:)); b = max(All(:));
xdata = 0:0.05:1;
zdata1 = All(1:21,1:21); zdata2 = All(1:21,22:42); zdata3 = All(1:21,43:63);
zdata4 = All(22:42,1:21); zdata5 = All(22:42,22:42); zdata6 = All(22:42,43:63);
zdata7 = All(43:63,1:21); zdata8 = All(43:63,22:42); zdata9 = All(43:63,43:63);
Figtype1(a,b,xdata,zdata1,zdata2,zdata3,zdata4,zdata5,zdata6,zdata7,zdata8,zdata9);