clearvars
close all
colormapfontsize=12;
axisfontsize=16;
panelletterfontsize=18;
figfolder= 'C:\Users\matlh\OneDrive - University of Leeds\My Projects\In-Leeds\Cooperator-Treshold system\Final Figure 4\';
panelletter={'(a)','(b)','(c)','(d)','(e)','(f)','(g)','(h)','(i)'};
delta=0;%0.32;
sS=.1; aAmM=.25; kKminus=120; kKplus=1000; meanK=(kKplus*(1+delta)+kKminus*(1-delta))/2;
cCurlyK=2*kKminus*kKplus./(kKplus*(1-delta)+kKminus*(1+delta));
figidx = {12,13,14};
N_TH=[60,80,100]; aux=cellfun(@(x)cat(2,num2str(x),'_'),num2cell(N_TH),'un',0); aux = [aux{:}]; aux = aux(1:end-1);
figname=['FinalFig4'];%_s' num2str(sS) '_a' num2str(aAmM) '_Kmn' num2str(kKminus) '_Kmx' num2str(kKplus) '_nth' aux];
loadfigdatafolder='C:\Users\matlh\OneDrive - University of Leeds\My Projects\In-Leeds\Cooperator-Treshold system\SimulationData\';
nu=[1e-4,1e-1,1e2];
% Initialize figure
f = figure('position',1e3*[0.1474    0.0490    1.0792    0.5008],'Name',figname);
set(f,'PaperPositionMode','auto');
set(f,'PaperOrientation','landscape');
ax0(1) = subplot(2,3,1,'parent',f); 
ax0(2) = subplot(2,3,2,'parent',f); 
ax0(3) = subplot(2,3,3,'parent',f); 
ax0(4) = subplot(2,3,4,'parent',f); 
ax0(5) = subplot(2,3,5,'parent',f); 
ax0(6) = subplot(2,3,6,'parent',f);
%Move axes. 
%ax0(5).Position=ax0(5).Position+[0.01 0 0 0];

%% Load and preprocess Matt's data for the "N_th_idx" N_th value.
N_th_idx=2;
nu_npoints=101; delta_npoints=nu_npoints;
data=nan(nu_npoints,delta_npoints,3,1e3); %dimensions: (N_TH,delta,nu,{N_C,N_D,t_fix},realization).
for j=1:1e3
    Aux = squeeze(load([loadfigdatafolder 'FixAndCoexProbData\heatmap K_0=560, gamma=0.79, Nth_values=[60, 80, 100], sims=1 ' num2str(j) '.mat'],'data'));
    data(:,:,:,j) = squeeze(Aux.data(1,N_th_idx,:,:,:)); %first dimension is (s,a) parameters, second and onwards as in "data" variable.
    clearvars Aux
end
%Set data in a proper format. Matt's simulations always ran up to t=2<K(nu=0,delta)> even with prior fixation.
N_th=N_TH(N_th_idx);
Data.nu = 10.^((-4)+(0:((nu_npoints)-1))'*((2)-(-4))/((nu_npoints)-1)); Data.nu=repmat(Data.nu,[delta_npoints 1]);
Data.delta = (-1)+(0:(delta_npoints-1))'*2/(delta_npoints-1); Data.delta=repelem(Data.delta,nu_npoints);
Data.N_C=cell(length(Data.nu),1); Data.N_D=Data.N_C; Data.t=Data.N_C;
for j=1:length(Data.nu)
    ij=1+floor((j-1)/delta_npoints); jj=1+rem(j-1,nu_npoints);
    Data.N_C{j} = squeeze(data(ij,jj,1,:));
    Data.N_D{j} = squeeze(data(ij,jj,2,:));
    Data.t{j} = squeeze(data(ij,jj,3,:));
end
clearvars data
%Process data.
Data.N=cellfun(@(x,y)x+y,Data.N_C,Data.N_D,'un',0);

%% Plot the first histogram in subplot m=1;
m=1; axes(ax0(m));
%Load and preprocess Matt's data for the distribution at nu=1e-4.
Y=readmatrix([loadfigdatafolder 'Distributions\distribution [N_th, nu, delta]=[80, ' num2str(nu(m)) ', 0.0].csv']); %load Matt's data.
N=sum(Y,2); %get N=N_C+N_D.
Sigma=10; %width of the Gaussian Process for the histogram.
x=0:round(kKplus+6*sqrt(kKplus),2,'significant'); %choose the N-range.
%Compute the N-histogram with a Gaussian Process.
n=zeros(size(x));
for i=1:numel(N)
n=n+exp(-(x-N(i)).^2/(2*Sigma^2))/(Sigma*sqrt(2*pi));
end
n=n/sum(n)*numel(x); %Normalize the result.
%Compute the N_C-histogram with a Gaussian Process.
nc=zeros(size(x));
for i=1:numel(N)
nc=nc+exp(-(x-Y(i,1)).^2/(2*Sigma^2))/(Sigma*sqrt(2*pi));
end
nc=nc/sum(nc)*numel(x); %Normalize the result.
%Compute the N_D-histogram with a Gaussian Process.
nd=zeros(size(x));
for i=1:numel(N)
nd=nd+exp(-(x-Y(i,2)).^2/(2*Sigma^2))/(Sigma*sqrt(2*pi));
end
nd=nd/sum(nd)*numel(x); %Normalize the result.
%Plot
plot(x,nc,'b-','LineWidth',2), hold on;
plot(x,nd,'r-','LineWidth',2), hold on;%plot([N_th-1,x+N_th],[0,nd],'r:','LineWidth',2), hold on;
plot(x,n,'-','Color',[0 .5 0],'LineWidth',2), hold on;
xlim([0 x(end)]); xlabel('# of individuals');
xticks([kKminus kKplus]); xticklabels(strsplit(num2str(xticks)));
ylabel('Pop. distr.'); yticklabels({});
ax0(m).Color='none';
set(gca,'FontSize',axisfontsize,'YScale','lin')
% Set the legend.
text(ax0(m),3/20,15/60,'N','Units','Normalized','HorizontalAlignment','left',...
        'VerticalAlignment','middle','Color',[0,.5,0],'FontWeight','Bold','FontSize',panelletterfontsize)
text(ax0(m),7/10,22/60,'N_{S}','Units','Normalized','HorizontalAlignment','left',...
        'VerticalAlignment','middle','Color',[1,0,0],'FontWeight','Bold','FontSize',panelletterfontsize)
text(ax0(m),1/20,35/60,'N_{R}','Units','Normalized','HorizontalAlignment','left',...
        'VerticalAlignment','middle','Color',[0,0,1],'FontWeight','Bold','FontSize',panelletterfontsize)
%Panel letter.
text(ax0(m),-.05,1.025,panelletter{m},'Units','Normalized','HorizontalAlignment','right',...
    'VerticalAlignment','bottom','Color',[0,0,0],'FontWeight','Bold','FontSize',panelletterfontsize)
title(ax0(m),['\nu = 10^{' num2str(log10(nu(m))) '}'],'FontSize',panelletterfontsize,'Units','Normalized','Position',[0.5 0.85])
box off
ax0(m).Position = ax0(m).Position + [0 0.06 0 0];

%% Plot the second histogram in subplot m=2;
m=2; axes(ax0(m));
%Load and preprocess Matt's data for the distribution at nu=1e-1.
Y=readmatrix([loadfigdatafolder 'Distributions\distribution [N_th, nu, delta]=[80, ' num2str(nu(m)) ', 0.0].csv']); %load Matt's data.
N=sum(Y,2); %get N=N_C+N_D.
Sigma=10; %width of the Gaussian Process for the histogram.
x=0:round(kKplus+6*sqrt(kKplus),2,'significant'); %choose the N-range.
%Compute the N-histogram with a Gaussian Process.
n=zeros(size(x));
for i=1:numel(N)
n=n+exp(-(x-N(i)).^2/(2*Sigma^2))/(Sigma*sqrt(2*pi));
end
n=n/sum(n)*numel(x); %Normalize the result.
%Compute the N_C-histogram with a Gaussian Process.
nc=zeros(size(x));
for i=1:numel(N)
nc=nc+exp(-(x-Y(i,1)).^2/(2*Sigma^2))/(Sigma*sqrt(2*pi));
end
nc=nc/sum(nc)*numel(x); %Normalize the result.
%Compute the N_D-histogram with a Gaussian Process.
nd=zeros(size(x));
for i=1:numel(N)
nd=nd+exp(-(x-Y(i,2)).^2/(2*Sigma^2))/(Sigma*sqrt(2*pi));
end
nd=nd/sum(nd)*numel(x); %Normalize the result.
%Plot
plot(x,nc,'b-','LineWidth',2), hold on;
plot(x,nd,'r-','LineWidth',2), hold on;%plot([N_th-1,x+N_th],[0,nd],'r:','LineWidth',2), hold on;
plot(x,n,'-','Color',[0 .5 0],'LineWidth',2), hold on;
xlim([0 x(end)]); xticks([kKminus kKplus]); xticklabels(strsplit(num2str(xticks))); yticklabels({});
ax0(m).Color='none';
set(gca,'FontSize',axisfontsize,'YScale','lin')
%Panel letter.
text(ax0(m),-.05,1.025,panelletter{m},'Units','Normalized','HorizontalAlignment','right',...
    'VerticalAlignment','bottom','Color',[0,0,0],'FontWeight','Bold','FontSize',panelletterfontsize)
title(ax0(m),['10^{' num2str(log10(nu(m))) '}'],'FontSize',panelletterfontsize,'Units','Normalized','Position',[0.5 0.85])
box off
ax0(m).Position = ax0(m).Position + [0 0.06 0 0];

%% Plot the third histogram in subplot m=3;
m=3; axes(ax0(m));
%Load and preprocess Matt's data for the distribution at nu=1e2.
Y=readmatrix([loadfigdatafolder 'Distributions\distribution [N_th, nu, delta]=[80, ' num2str(nu(m)) ', 0.0].csv']); %load Matt's data.
N=sum(Y,2); %get N=N_C+N_D.
Sigma=10; %width of the Gaussian Process for the histogram.
x=0:round(kKplus+6*sqrt(kKplus),2,'significant'); %choose the N-range.
%Compute the N-histogram with a Gaussian Process.
n=zeros(size(x));
for i=1:numel(N)
n=n+exp(-(x-N(i)).^2/(2*Sigma^2))/(Sigma*sqrt(2*pi));
end
n=n/sum(n)*numel(x); %Normalize the result.
%Compute the N_C-histogram with a Gaussian Process.
nc=zeros(size(x));
for i=1:numel(N)
nc=nc+exp(-(x-Y(i,1)).^2/(2*Sigma^2))/(Sigma*sqrt(2*pi));
end
nc=nc/sum(nc)*numel(x); %Normalize the result.
%Compute the N_D-histogram with a Gaussian Process.
nd=zeros(size(x));
for i=1:numel(N)
nd=nd+exp(-(x-Y(i,2)).^2/(2*Sigma^2))/(Sigma*sqrt(2*pi));
end
nd=nd/sum(nd)*numel(x); %Normalize the result.
%Plot
plot(x,nc,'b-','LineWidth',2), hold on;
plot(x,nd,'r-','LineWidth',2), hold on;%plot([N_th-1,x+N_th],[0,nd],'r:','LineWidth',2), hold on;
plot(x,n,'-','Color',[0 .5 0],'LineWidth',2), hold on;
xlim([0 2*round(cCurlyK)]); xticks([0 round(cCurlyK)]); xticklabels(strsplit(num2str(xticks))); yticklabels({}); ylim([0 2*max(n)]);
ax0(m).Color='none';
set(gca,'FontSize',axisfontsize,'YScale','lin')
text(ax0(m),1.1*round(cCurlyK)/max(xlim),41/60,'N = ','Units','Normalized','HorizontalAlignment','right',...
        'VerticalAlignment','middle','Color',[0,.5,0],'FontWeight','Bold','FontSize',panelletterfontsize)
text(ax0(m),1.1*round(cCurlyK)/max(xlim),40/60,'$\mathcal{K}$','Units','Normalized','HorizontalAlignment','left',...
        'VerticalAlignment','middle','Color',[0,.5,0],'FontWeight','Bold','FontSize',panelletterfontsize,'Interpreter','latex')
%Panel letter.
text(ax0(m),-.05,1.025,panelletter{m},'Units','Normalized','HorizontalAlignment','right',...
    'VerticalAlignment','bottom','Color',[0,0,0],'FontWeight','Bold','FontSize',panelletterfontsize)
title(ax0(m),['10^{' num2str(log10(nu(m))) '}'],'FontSize',panelletterfontsize,'Units','Normalized','Position',[0.5 0.85])
box off
ax0(m).Position = ax0(m).Position + [0 0.06 0 0];

%% Plot the <N> vs nu in subplot m=4;
m=4; axes(ax0(m));
%Compute the data for this plot.
idx=abs(Data.delta-delta)<1e-6;
x=Data.nu(idx);
y(:,1)=cellfun(@mean,Data.N_C(idx),'un',1);
y(:,2)=cellfun(@mean,Data.N_D(idx),'un',1);
y(:,3)=cellfun(@mean,Data.N(idx),'un',1);
%Gaussian-filter data.
y_smooth=10.^(smoothdata(log10(y),1,'gaussian',10));
%Plot.
plot(x,y_smooth(:,1),'b-','LineWidth',2), hold on;
plot(x,y_smooth(:,2),'r-','LineWidth',2);
plot(x,y_smooth(:,3),'-','Color',[0 .5 0],'LineWidth',2);
%plot(x,y(:,1),'b:','LineWidth',2), hold on;
%plot(x,y(:,2),'r:','LineWidth',2);
%plot(x,y(:,3),'-','Color',[0 .5 0],'LineWidth',2);
plot([1e-4 1e-1],[meanK meanK],':','Color',[0 .5 0],'LineWidth',1.5);
%plot([1e0 1e2],[cCurlyK cCurlyK],'k:','LineWidth',1.5);
xlim([1e-4 1e2]); xticks(10.^(-4:3:2)); xlabel('\nu (a.u.)');
yticks([0 round(cCurlyK) round(meanK)]);
ax0(m).Color='none';
set(gca,'FontSize',axisfontsize,'XScale','log')
%Legend and additions
text(ax0(m),0.5,1,'$<K>$','Units','Normalized','HorizontalAlignment','left',...
        'VerticalAlignment','top','Color',[0,.5,0],'FontWeight','Bold','FontSize',14/18*panelletterfontsize,'Interpreter','latex')
text(ax0(m),0.6,1.1*mean([meanK,cCurlyK])/max(ylim),'<N>','Units','Normalized','HorizontalAlignment','left',...
        'VerticalAlignment','middle','Color',[0,.5,0],'FontWeight','Bold','FontSize',panelletterfontsize)
text(ax0(m),0.90,1.25*cCurlyK/max(ylim),'$\mathcal{K}$','Units','Normalized','HorizontalAlignment','left',...
        'VerticalAlignment','middle','Color',[0,.5,0],'FontWeight','Bold','FontSize',14/18*panelletterfontsize,'Interpreter','latex')
text(ax0(m),0.325,1.1*mean([meanK,cCurlyK])/max(ylim),'<N_R>','Units','Normalized','HorizontalAlignment','center',...
        'VerticalAlignment','middle','Color','b','FontWeight','Bold','FontSize',panelletterfontsize)
text(ax0(m),0.3,0.3*mean([meanK,cCurlyK])/max(ylim),'<N_S>','Units','Normalized','HorizontalAlignment','center',...
        'VerticalAlignment','middle','Color','r','FontWeight','Bold','FontSize',panelletterfontsize)
text(ax0(m),1.12,1.25*(cCurlyK-N_th)/max(ylim),'$\mathcal{K}-N_{th}$','Units','Normalized','HorizontalAlignment','right',...
        'VerticalAlignment','middle','Color','r','FontWeight','Bold','FontSize',14/18*panelletterfontsize,'Interpreter','latex')
text(ax0(m),0.99,0.55*N_th/max(ylim),'$N_{th}$','Units','Normalized','HorizontalAlignment','right',...
        'VerticalAlignment','middle','Color','b','FontWeight','Bold','FontSize',14/18*panelletterfontsize,'Interpreter','latex')
%Panel letter.
text(ax0(m),-.05,1.025,panelletter{m},'Units','Normalized','HorizontalAlignment','right',...
    'VerticalAlignment','bottom','Color',[0,0,0],'FontWeight','Bold','FontSize',panelletterfontsize)
box off
ax0(m).Position = ax0(m).Position.*[1 1 0.8 1] + [0 0 0 0];

%% Plot the phase diagram in subplot m=5;
m=5; axes(ax0(m));
%Plot data.
Color=cell2mat(cellfun(@(x)[0,mean(x),0],Data.N,'un',0)); Color=Color/max(max(Color));
scatter(ax0(m),Data.nu,Data.delta,20,Color,'s','filled'); set(gca,'XScale','log');
xlim(10.^[-4,2].*10.^[-.5,.5])
ylim([-1 1]+[-1 1]/6)
xticks(10.^(-4:3:2))
%xlabel('\nu (a.u.)')
yticks(-1:1)
ylabel('\delta')
title(['<N>'])
ax0(m).Color='none';
set(gca,'FontSize',axisfontsize)
axis square
colormap(ax0(m),sort(Color));
aux = colorbar(ax0(m),'TickLabels',{'0','5','10'}); aux.Title.String='    x10^2'; aux.Title.FontSize=aux.Title.FontSize*9/10;
set(gca,'Fontsize',axisfontsize)
text(ax0(m),-.05...-0.15,
    ,1.025,panelletter{m},'Units','Normalized','HorizontalAlignment','right',...
    'VerticalAlignment','bottom','Color',[0,0,0],'FontWeight','Bold','FontSize',panelletterfontsize)
ax0(m).Position = ax0(m).Position + [-0.06 0 0 0];

%% Plot the phase diagram in subplot m=6;
m=6; axes(ax0(m));
%Plot data.
CoopFrac=cellfun(@(x,y)x./(x+y),Data.N_C,Data.N_D,'un',0);
NoCoexistence=cellfun(@(x,y)x<min(2*mean(y),max(x)),Data.t,Data.N,'un',0); %cellfun(@(x,y,z)x<min(2*mean(y),z),Data.t,Data.N,num2cell(maxt),'un',0);
Color=cell2mat(cellfun(@(x,y)mean(y.*[1,1,1]+(1-y).*[1-x,0*x,x]),CoopFrac,NoCoexistence,'un',0));
scatter(ax0(m),Data.nu,Data.delta,20,Color,'s','filled'); set(gca,'XScale','log');
xlim(10.^[-4,2].*10.^[-.5,.5])
ylim([-1 1]+[-1 1]/6)
xticks(10.^(-4:3:2))
%xlabel('\nu (a.u.)')
yticks(-1:1)
%ylabel('\delta')
%title(['<x>'])
ax0(m).Color='none';
set(ax0(m),'FontSize',axisfontsize)
axis square
ax0(m).Position = ax0(m).Position + [-0.08 0 0 0];
%Create ternary plot with colormap.
auxax=axes(f,'Position',[0.84 0.165 0.075 0.225],'Units','Normalized');
Vertex=[0,0;0,1;sin(pi/3),0.5];
VertexColor=[[0,0,1];[1,0,0];[1,1,1]];
Npoints=200+1; aux=((0:(Npoints-1))/(Npoints-1))'.*[1,1];
Points(:,1)=repmat(aux(:,1),[Npoints 1]); Points(:,2)=repelem(aux(:,2),Npoints);
Slope = (Vertex(3,2)-Vertex(1,2))/(Vertex(3,1)-Vertex(1,1));
Points=Points((Points(:,1)*Slope<=Points(:,2))&((1-Points(:,1)*Slope)>=Points(:,2)),:);
Color=nan(size(Points,1),3);
for i=1:size(Points,1)
    Color(i,:)=[1-Points(i,2)-Points(i,1)*Slope,Points(i,2)-Points(i,1)*Slope,Points(i,1)/Vertex(3,1)]; Color(i,:)=sum(Color(i,:)'.*VertexColor,1);
end
scatter(auxax,Points(:,1),Points(:,2),10,Color,'>','filled'); axis square;
set(gca,'Color','none','XColor','none','YColor','none');
text(gca,0,0,{'Resist.';'coex. frac.'},'Units','Normalized','HorizontalAlignment','center',...
    'VerticalAlignment','top','Color',VertexColor(1,:),'FontWeight','Bold','FontSize',colormapfontsize,'Rotation',180/pi*atan(-Slope))
text(gca,0,1,{'Sens.';'coex. frac.'},'Units','Normalized','HorizontalAlignment','center',...
    'VerticalAlignment','bottom','Color',VertexColor(2,:),'FontWeight','Bold','FontSize',colormapfontsize,'Rotation',180/pi*atan(Slope))
text(gca,sin(pi/3),0.5,{'Any fix.'},'Units','Normalized','HorizontalAlignment','center',...
    'VerticalAlignment','bottom','Color',[0,0,0],'FontWeight','Bold','FontSize',colormapfontsize,'Rotation',-90)
text(ax0(m),-.05,1.025,panelletter{m},'Units','Normalized','HorizontalAlignment','right',...
    'VerticalAlignment','bottom','Color',[0,0,0],'FontWeight','Bold','FontSize',panelletterfontsize)
auxax.Position = auxax.Position + [0 0 0 0];
% Put the reference of the environments of d-f in panels a-c.
panelletter2 = {'  (a)','(b)','(c)  '}; HorzAlign = {'left','center','right'};
textposition(:,:,1) = [1,0;1,0.15;1,0]; textposition(:,:,2) = [1,0;1,0.15;1,0];
Clors(:,:,1) = [1,1,1;1,1,1;1,1,1]; Clors(:,:,2) = [0,0,0;0,0,0;0,0,0];

for i=1:3
    axes(ax0(5)); hold on; scatter(ax0(5),nu(i),0,100,Clors(i,:,1),'Marker','*','Linewidth',2);
    text(ax0(5),nu(i)*textposition(i,1,1),0+textposition(i,2,1),panelletter2{i},'HorizontalAlignment',HorzAlign{i},...
        'VerticalAlignment','bottom','Color',Clors(i,:,1),'FontWeight','Bold','FontSize',axisfontsize)
    axes(ax0(6)); hold on; scatter(ax0(6),nu(i),0,100,Clors(i,:,2),'Marker','*','Linewidth',2);
    text(ax0(6),nu(i)*textposition(i,1,2),0+textposition(i,2,2),panelletter2{i},'HorizontalAlignment',HorzAlign{i},...
        'VerticalAlignment','bottom','Color',Clors(i,:,2),'FontWeight','Bold','FontSize',axisfontsize)    
end

%% Save Figure
set(gcf, 'InvertHardCopy', 'off'); set(gcf,'Color','w');  %So that the (a), (b), (c) text is not in an inverted color with respect to the background when saving the image.
savefig(f,[figfolder figname '.fig'])
saveas(f,[figfolder figname '.png'],'png')
saveas(f,[figfolder figname '.pdf'],'pdf')
