clearvars
close all
colormapfontsize=14;
axisfontsize=16;
panelletterfontsize=18;
figfolder= 'C:\Users\matlh\OneDrive - University of Leeds\My Projects\In-Leeds\Cooperator-Treshold system\Final Figure 2\';
panelletter={'(a)','(b)','(c)','(d)','(e)','(f)'};
sS=.1; aAmM=.25; kKminus=120; kKplus=1000;
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=['FinalFig2'];
loadfigdatafolder='C:\Users\matlh\OneDrive - University of Leeds\My Projects\In-Leeds\Cooperator-Treshold system\SimulationData\';
Traj=[9,5,1];
% Initialize figure
f = figure('position',[13 53 1325.1 595.2],'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);
aux_x = 0.05; aux_y = .15;
ax0(1).Position=ax0(1).Position+[-0*aux_x -aux_y -aux_x +aux_y];
ax0(2).Position=ax0(2).Position+[-1*aux_x -aux_y -aux_x +aux_y];
ax0(3).Position=ax0(3).Position+[-2*aux_x -aux_y -aux_x +aux_y];
ax0(4).Position=ax0(4).Position+[0 0 0 -aux_y];
ax0(5).Position=ax0(5).Position+[0 0 0 -aux_y];
ax0(6).Position=ax0(6).Position+[0 0 0 -aux_y];
titl(1:3)={''}; titl{1}='N_{th} = ';
% Plot the phase diagrams
m=0;
for i=[figidx{:}]%11:14
    m=m+1;
    %Load Matt's data for the "m" N_th value.
    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,m,:,:,:)); %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(m);
    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.
    maxt=2*(kKminus*(1-Data.delta)/2+kKplus*(1+Data.delta)/2);
    Data.N=cellfun(@(x,y)x+y,Data.N_C,Data.N_D,'un',0);
    NoCoexistence=cellfun(@(x,y,z)x<min(2*mean(y),max(x)),Data.t,Data.N,num2cell(maxt),'un',0); %cellfun(@(x,y,z)x<min(2*mean(y),z),Data.t,Data.N,num2cell(maxt),'un',0);
    CoopFrac=cellfun(@(x,y)x./(x+y),Data.N_C,Data.N_D,'un',0);
    Color=cell2mat(cellfun(@(x,y)[mean((1-x).*y),0,mean(x.*y)],CoopFrac,NoCoexistence,'un',0));
    Color2=cell2mat(cellfun(@(x)[0,mean(x),0],Data.N,'un',0)); Color2=Color2/max(max(Color2));
    g(m) = figure; set(gcf,'position',[500   100   1.25*560   1.25*420],'Name',['AsymDensDepSys' num2str(i) '_CoexAndDom'])
    scatter(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:2:2))
    xlabel('\nu (a.u.)')
    yticks(-1:.5:1)
    ylabel('\delta')
    title(['s=',num2str(sS),', a=',num2str(aAmM),', N_{Th}=',num2str(N_TH(m)),', K^-=',num2str(kKminus),', K^+=',num2str(kKplus)])
    set(gca,'FontSize',axisfontsize)
    axis square
    %Copying objects to final figure.
    g(m)=gcf;
    ax = gca;    
    axcp(m) = copyobj(ax, f);
    set(axcp(m),'Position',get(ax0(m),'position'));
    delete(ax0(m));    
    sct=findobj(axcp(m),'Type','Scatter');
    sct.SizeData=75;
    figure(f)
    s=sS; a=aAmM; N_th=N_TH(m);
    %plot the contour for the optimal environments to eradicate resistant cooperators with transient cooperator dips.
    set(f,'CurrentAxes',axcp(m)); hold on;
    tdip=log((1-s)/(a-s)*(1-kKminus/kKplus));
    nu=10.^(-4:(2-(-4))/2400:2); alph=2*log((kKplus)/(kKminus));%=4;
    y(1,:)=1-s./(alph*nu); y(2,:)=1./(tdip*nu)-1; y(3,:)=1-s./(nu); y(4,:)=(a-s)./(alph*nu)-1;
    yy{1}=min(y(1:2,:)); yy{2}=max(y(3:4,:)); clearvars y
    y(1,:)=yy{1}(yy{1}>=yy{2}); y(2,:)=yy{2}(yy{1}>=yy{2}); nu=nu(yy{1}>=yy{2});
    y=cat(2,y(1,:),flip(y(2,:)),y(1,1)); nu=cat(2,nu,flip(nu),nu(1));
    plot(axcp(m),nu,y,'g-','LineWidth',2); clearvars y
    axcp(m).Title.String = [titl{m},num2str(N_th)]; yticks(-1:1); xticks([1e-4 1e-1 1e2]); xlabel('\nu (a.u.)');
    if m>1; ylabel(''); xlabel(''); end%xticklabels({}); end
    set(gca,'Fontsize',axisfontsize)
    text(axcp(m),-.05,1.05,panelletter{m},'Units','Normalized','HorizontalAlignment','right',...
        'VerticalAlignment','bottom','Color',[0,0,0],'FontWeight','Bold','FontSize',panelletterfontsize)
end

%Create ternary plot with colormap.
auxax=axes(f,'Position',[0.825 0.55 0.075 0.225],'Units','Normalized');
Vertex=[0,0;0,1;sin(pi/3),0.5];
VertexColor=[[0,0,1];[1,0,0];[0,0,0]];
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.';'fix.'},'Units','Normalized','HorizontalAlignment','center',...
    'VerticalAlignment','top','Color',VertexColor(1,:),'FontWeight','Bold','FontSize',colormapfontsize,'Rotation',180/pi*atan(-Slope))
text(gca,0,1,{'Sens.';'fix.'},'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,{'Coex.'},'Units','Normalized','HorizontalAlignment','center',...
    'VerticalAlignment','bottom','Color',VertexColor(3,:),'FontWeight','Bold','FontSize',colormapfontsize,'Rotation',-90)
%title(gca,{'Probability of'},'Units','Normalized','Position',[Vertex(3,1)/2 1.1],'FontWeight','Bold','FontSize',panelletterfontsize)

% Plot the trajectories
traj(m+(1:3))=Traj;
N_Th(m+(1:3))=N_TH([2 2 2]);
nu_delta(m+1,:)=[1e-3,-.5]; nu_delta(m+2,:)=[1e-1,2/3]; nu_delta(m+3,:)=[2e1,1/3]; % example environments
AuxDeltaStr{m+1,1}=num2str(nu_delta(m+1,1)); AuxDeltaStr{m+2,1}=num2str(nu_delta(m+2,1)); AuxDeltaStr{m+3,1}='20.0';
AuxDeltaStr{m+1,2}=num2str(nu_delta(m+1,2)); AuxDeltaStr{m+2,2}='0.6666666666666666'; AuxDeltaStr{m+3,2}='0.3333333333333333';
mm=0;
for i=[figidx{:}]%11:14
    m=m+1;
    mm=mm+1;
    %load data
    %load Matt's data and set it in a proper format.
    TrajData=cell(10,1); %Here we set the number of realizations (10).
    N_th=N_Th(m);
    ll=0;
    for j=0:(size(TrajData,1)-1)
        ll=ll+1;
        TrajData{ll,1}(:,:)=readmatrix([loadfigdatafolder 'Trajectories\trajectory [N_th, nu, delta]=[' num2str(N_th) ', ' AuxDeltaStr{m,1} ', ' AuxDeltaStr{m,2} '] ' num2str(j) '.csv'],'NumHeaderLines',1);
    end
    TrajData_Size = cellfun(@(x)size(x,1),TrajData,'un',1); Max_TrajData_Size = max(TrajData_Size);
    t_all = {nan(Max_TrajData_Size,ll)}; k_all = t_all; Nc_all = t_all; Nd_all = t_all; %Initialize variables.
    for j=1:ll %One-element cell with a reactiontimes-by-realization matrix in it, filled with nans when a realization finished first, where the values correspond to the name of the variable "*_all".
        Nc_all{1}(1:TrajData_Size(j),j) = TrajData{j}(:,1);
        Nd_all{1}(1:TrajData_Size(j),j) = TrajData{j}(:,2);
        k_all{1}(1:TrajData_Size(j),j) = TrajData{j}(:,3);
        t_all{1}(1:TrajData_Size(j),j) = TrajData{j}(:,4);
    end
    clearvars TrajData TrajData_Size Max_TrajData_Size ll j
    Pixel=1;
    %plot
    for ij=Pixel % in case we want to explore more pixels.
        g(m) = figure; set(gcf,'position',[33.8 67.4 1439.2 703.2],'Name',['AsymDensDepSys' num2str(i) '_Traj_pixel' num2str(ij)])
        for jj=1:size(t_all{ij},2)
            Time = t_all{ij}(:,jj); NC = Nc_all{ij}(:,jj); Aux = Nd_all{ij}(:,jj);
            ND1 = nan(size(Aux)); ND2 = ND1; ND1(NC<N_th) = Aux(NC<N_th); ND2(NC>=N_th) = Aux(NC>=N_th);
            Kval = k_all{ij}(:,jj);
            nsubrow=3; nsubcol=4;
            ax1(jj)=subplot(nsubrow,nsubcol,jj); hold on;
            plot(Time,N_th*ones(size(Time)),'b--')
            plot(Time,Kval,'k')
            plot(Time,NC,'b')
            plot(Time,ND1,'r')
            plot(Time,ND2,'g')
        end
    end
    %traj(m)=2;
    ax=ax1(traj(m));
    axcp(m) = copyobj(ax, f);
    set(axcp(m),'Position',get(ax0(m),'position'));
    delete(ax0(m));
    figure(f)
    set(f,'CurrentAxes',axcp(m));
    % One single red line for simulated Nc.
    o=findobj(gca,'Type','Line',{'Color',[0,1,0],'-or','Color',[1,0,0]});
    o(1).YData=nanmean(cat(1,reshape(o(1).YData,1,[]),reshape(o(2).YData,1,[])));
    set(o(1),'Color',[1,0,0]); delete(o(2));
    ND=cat(1,reshape(o(1).XData,1,[]),reshape(o(1).YData,1,[]));
    o=findobj(gca,'Type','Line','Color',[0,0,1],'LineStyle','-');
    NC=cat(1,reshape(o(1).XData,1,[]),reshape(o(1).YData,1,[]));
    % Thicken deterministic lines (K).
    gto_aux=findobj(gca,'Type','Line','Color',[0,0,0]);
    xlim_aux=500;
    if m<6
        set(gto_aux,'LineWidth',2);
    else
        XdAta=get(gto_aux,'XData'); YdAta=get(gto_aux,'YData'); delta_aux=nu_delta(m,2); Kminus_aux=min(YdAta); Kplus_aux=max(YdAta);
        CurlyK_aux=2*Kminus_aux*Kplus_aux./(Kplus_aux*(1-delta_aux)+Kminus_aux*(1+delta_aux));
        set(gto_aux,'LineWidth',1e-6,'LineStyle',':','XData',[XdAta(1),xlim_aux],'YData',CurlyK_aux*[1,1],'LineWidth',2);
        set(gca,'YLim',[min(ylim),CurlyK_aux*1.5/1.025]);
    end
    % Aesthetics
    axcp(m).Title.String = ['\nu = ' num2str(nu_delta(m,1)) ' \delta = ' num2str(nu_delta(m,2))];
    xlabel('Time (a.u.)'); ylabel('Individuals');
    xLim=xlim; xlim([xLim(1) min(xLim(2),xlim_aux)]);
    if m<6; yticks([0 kKminus kKplus]); else; yticks([0 kKminus round(CurlyK_aux)]); end
     yticklabels(strsplit(num2str(yticks)));
    if m>4; ylabel(''); xlabel(''); end
    set(gca,'Fontsize',axisfontsize);
    set(gca,'YLim',1.025*get(gca,'YLim'))
    axcp(m).Title.Position=axcp(m).Title.Position.*[1 .85 1];
    axcp(m).Title.Visible='off';
    text(axcp(m),-.05,1.05,panelletter{m},'Units','Normalized','HorizontalAlignment','right',...
        'VerticalAlignment','bottom','Color',[0,0,0],'FontWeight','Bold','FontSize',panelletterfontsize)
    % clear big variables
    clearvars t_all k_all Nc_all Nd_all K
end
% Put the reference of the environments of d-f in panels a-c.
panelletter2={'(d) ','(e)  ','(f)  '}; HorzAlign = {'right','right','right'};
for i=1:3
    hold on; scatter(axcp(2),nu_delta(i+3,1),nu_delta(i+3,2),100,[1,1,1],'Marker','*','Linewidth',2);
    text(axcp(2),nu_delta(i+3,1),nu_delta(i+3,2),panelletter2{i},'HorizontalAlignment',HorzAlign{i},...
        'VerticalAlignment','bottom','Color',[1,1,1],'FontWeight','Bold','FontSize',axisfontsize)
end
% Set the legend for panel d and thus e-f.
text(axcp(4),0.2,50/60,'K','Units','Normalized','HorizontalAlignment','center',...
        'VerticalAlignment','bottom','Color',[0,0,0],'FontWeight','Bold','FontSize',panelletterfontsize)
text(axcp(4),0.6,40/max(axcp(4).YLim),'N_{S}','Units','Normalized','HorizontalAlignment','left',...
        'VerticalAlignment','middle','Color',[1,0,0],'FontWeight','Bold','FontSize',panelletterfontsize)
text(axcp(4),0.75,55/60,'N_{R}','Units','Normalized','HorizontalAlignment','center',...
        'VerticalAlignment','bottom','Color',[0,0,1],'FontWeight','Bold','FontSize',panelletterfontsize)
text(axcp(6),1.01,290/max(axcp(6).YLim),'$\mathcal{K}$','Units','Normalized','HorizontalAlignment','left',...
        'VerticalAlignment','middle','Color',[0,0,0],'FontWeight','Bold','FontSize',panelletterfontsize,'Interpreter','latex')

%% Save Figure
savefig(f,[figfolder figname '.fig'])
saveas(f,[figfolder figname '.png'],'png')
saveas(f,[figfolder figname '.pdf'],'pdf')
