clearvars
close all
colormapfontsize=14;
axisfontsize=16;
panelletterfontsize=18;
figfolder= 'C:\Users\matlh\OneDrive - University of Leeds\Research\My Projects\In-Leeds\Coop-FixRelThr\Research Data Leeds Repository\';
panelletter={'(a)','(b)','(c)','(d)','(e)','(f)','(g)','(h)','(i)','(j)','(k)','(l)'};
sS=.1; aAmM=.25; kKminus=100; kKplus=1000;
row_col = [3,4];
X_TH=[3,11,19,27]/30; aux=cellfun(@(x)cat(2,num2str(x),'_'),num2cell(X_TH),'un',0); aux = [aux{:}]; aux = aux(1:end-1);
figname=['FinalFig3'];
loadfigdatafolder='C:\Users\matlh\OneDrive - University of Leeds\Research\My Projects\In-Leeds\Coop-FixRelThr\Research Data Leeds Repository\SimulationData\';
Traj=[5,1,1,1];%[9,5,1,1];
% Initialize figure
f = figure('position',[13 53 1510 720],'Name',figname);%figure('position',[13 53 1325.1 595.2],'Name',figname);
set(f,'PaperPositionMode','auto');
set(f,'PaperOrientation','landscape');
ax0(1) = subplot(3,4,1,'parent',f); 
ax0(2) = subplot(3,4,2,'parent',f); 
ax0(3) = subplot(3,4,3,'parent',f); 
ax0(4) = subplot(3,4,4,'parent',f); 
ax0(5) = subplot(3,4,5,'parent',f); 
ax0(6) = subplot(3,4,6,'parent',f);
ax0(7) = subplot(3,4,7,'parent',f); 
ax0(8) = subplot(3,4,8,'parent',f); 
ax0(9) = subplot(3,4,9,'parent',f); 
ax0(10) = subplot(3,4,10,'parent',f); 
ax0(11) = subplot(3,4,11,'parent',f); 
ax0(12) = subplot(3,4,12,'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}='\{x_{th},\nu,\delta\} = ';

m=0;
traj(m+(1:4))=Traj;
X_Th(m+(1:4))=X_TH([1 2 3 4]);
AuxXthStr{m+1}= num2str(X_Th(m+1)); nu_delta(m+1,:)=[1e-4,-.5]; AuxDeltaStr{m+1,1}=num2str(nu_delta(m+1,1)); AuxDeltaStr{m+1,2}=num2str(nu_delta(m+1,2));
AuxXthStr{m+2}= ['0.36666666666666664']; nu_delta(m+2,:)=[1e-1,.5]; AuxDeltaStr{m+2,1}=num2str(nu_delta(m+2,1)); AuxDeltaStr{m+2,2}=num2str(nu_delta(m+2,2));
AuxXthStr{m+3}= ['0.6333333333333333']; nu_delta(m+3,:)=[1e-1,-.5]; AuxDeltaStr{m+3,1}=num2str(nu_delta(m+3,1)); AuxDeltaStr{m+3,2}=[num2str(nu_delta(m+3,2)) ''];
AuxXthStr{m+4}= num2str(X_Th(m+4)); nu_delta(m+4,:)=[1e2,.5]; AuxDeltaStr{m+4,1}=[num2str(nu_delta(m+4,1)) '.0']; AuxDeltaStr{m+4,2}=num2str(nu_delta(m+4,2)); % example environments

% Plot the first row: composition
mm=0;
for i=1:row_col(2)%11:14
    m=m+1;
    M=rem(m,row_col(2));  M(M==0)=4;
    mm=mm+1;
    %load data
    %load Matt's data and set it in a proper format.
    TrajData=cell(1,1); %cell(10,1); %Here we set the number of realizations (10).
    X_th=X_Th(M);
    ll=0;
    for j=0:(size(TrajData,1)-1)
        ll=ll+1;
        TrajData{ll,1}(:,:)=readmatrix([loadfigdatafolder 'Trajectories\trajectory [x_th, nu, delta]=[' AuxXthStr{M} ', ' AuxDeltaStr{M,1} ', ' AuxDeltaStr{M,2} '] ' num2str(traj(M)-1) '.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<X_th) = Aux(NC<X_th); ND2(NC>=X_th) = Aux(NC>=X_th);
            ND = nanmean(cat(2,ND1,ND2),2);
            Kval = k_all{ij}(:,jj);
            nsubrow=3; nsubcol=4;
            ax1(jj)=subplot(nsubrow,nsubcol,jj); hold on;
            if M==2||M==3 % plot the shaded background.
                aux_idx = find(~isnan(Kval),1,'last');
                if aux_idx~=size(Kval,1); Kval(aux_idx:end)=Kval(aux_idx); end                
                aux=logical(cat(1,1,abs(diff(Kval))/(max(Kval)-min(Kval)))+cat(1,abs(diff(Kval))/(max(Kval)-min(Kval)),1)); % logical vector of indeces prior and posterior to a K switch.
                Kval_aux = Kval(aux); Kval_aux = (Kval_aux-min(Kval_aux))/(max(Kval_aux)-min(Kval_aux)); Time_aux = Time(aux);
                %fill(cat(1,Time_aux,flip(Time_aux)),cat(1,Kval_aux,flip(0*Kval_aux)), 'k', 'edgecolor', 'none', 'facealpha', 0.2);
                area(Time_aux,Kval_aux,0,'LineStyle','none','FaceColor', 'k',  'facealpha', 0.2);
            end
            plot(Time,X_th*ones(size(Time)),'k--');
            plot(Time,NC./(NC+ND),'b')
        end
    end
    %traj(m)=2;
    ax=ax1(1);
    axcp(m) = copyobj(ax, f);
    set(axcp(m),'Position',get(ax0(m),'position'));
    delete(ax0(m));
    figure(f)
    set(f,'CurrentAxes',axcp(m));
    % Thicken deterministic lines (K).
    gto_aux=findobj(gca,'Type','Line','Color',[0,0,0]);
    xlim_aux=200;
    set(gto_aux,'LineWidth',2);
    % Aesthetics
    axcp(m).Title.String = ['\{' num2str(round(X_Th(M),2)) ', ' num2str(nu_delta(M,1)) ', ' num2str(nu_delta(M,2)) '\}'];
    %xlabel('Time (a.u.)');    
    ylim([0 1]); yTicks=[0 X_th 1]; yticks(yTicks); yticklabels(cellfun(@(x)num2str(round(x,2)),num2cell(yTicks),'un',0))   
    xLim=xlim; xlim([xLim(1) min(xLim(2),xlim_aux)]);
    if M==1
        axcp(m).Title.String = ['\{' num2str(round(X_Th(M),2)) ', ' '10^{-4}' ', ' num2str(nu_delta(M,2)) '\}'];
        axcp(m).Title.String = ['\{x_{th},\nu,\delta\} = ' axcp(m).Title.String];
        ylabel('Pop. comp. x');
    end
    if M>1; ylabel(''); axcp(m).Title.Position=axcp(m).Title.Position.*[1 1.05 1]; end
    set(gca,'Fontsize',axisfontsize);
    axcp(m).Title.FontWeight = 'Normal'; axcp(m).Title.FontSize = axisfontsize;
    axcp(m).Title.Visible='on';
    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

% Plot the second row: total population
mm=0;
for i=1:row_col(2)%11:14
    m=m+1;
    M=rem(m,row_col(2)); M(M==0)=row_col(2);
    mm=mm+1;
    %load data
    %load Matt's data and set it in a proper format.
    TrajData=cell(1,1); %cell(10,1); %Here we set the number of realizations (10).
    X_th=X_Th(M);
    ll=0;
    for j=0:(size(TrajData,1)-1)
        ll=ll+1;
        TrajData{ll,1}(:,:)=readmatrix([loadfigdatafolder 'Trajectories\trajectory [x_th, nu, delta]=[' AuxXthStr{M} ', ' AuxDeltaStr{M,1} ', ' AuxDeltaStr{M,2} '] ' num2str(traj(M)-1) '.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<X_th) = Aux(NC<X_th); ND2(NC>=X_th) = Aux(NC>=X_th);
            ND = nanmean(cat(2,ND1,ND2),2);
            Kval = k_all{ij}(:,jj);
            nsubrow=3; nsubcol=4;
            ax1(jj)=subplot(nsubrow,nsubcol,jj); hold on;
            plot(Time,Kval,'k')
            plot(Time,NC+ND,'Color',[0 .75 0])
        end
    end
    %traj(m)=2;
    ax=ax1(1);
    axcp(m) = copyobj(ax, f);
    set(axcp(m),'Position',get(ax0(m),'position'));
    delete(ax0(m));
    figure(f)
    set(f,'CurrentAxes',axcp(m));
    % Thicken deterministic lines (K).
    gto_aux=findobj(gca,'Type','Line','Color',[0,0,0]);
    xlim_aux=200;
    if nu_delta(M,1)<1e2
        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,'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 = ['\{' num2str(round(X_Th(M),2)) ', ' num2str(nu_delta(M,1)) ', ' num2str(nu_delta(M,2)) '\}'];
    %xlabel('Time (a.u.)');
    if M==1
        axcp(m).Title.String = ['\{' num2str(round(X_Th(M),2)) ', ' '10^{-4}' ', ' num2str(nu_delta(M,2)) '\}'];
        axcp(m).Title.String = ['\{x_{th},\nu,\delta\} = ' axcp(m).Title.String];
        ylabel('Total pop. N'); ylim([0 kKminus*2/1.25]);
    end
    xLim=xlim; xlim([xLim(1) min(xLim(2),xlim_aux)]);
    if M~=row_col(2); yticks([0 kKminus kKplus]); else; yticks([0 kKminus round(CurlyK_aux)]); end
    yticklabels(strsplit(num2str(yticks)));
    if M>1; ylabel(''); axcp(m).Title.Position=axcp(m).Title.Position.*[1 1.1 1]; end
    set(gca,'Fontsize',axisfontsize);
    set(gca,'YLim',1.025*get(gca,'YLim')); ylim([0 max(ylim)]);
    axcp(m).Title.FontWeight = 'Normal'; axcp(m).Title.FontSize = axisfontsize;
    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

% Plot the third row: individuals
mm=0;
for i=1:row_col(2)%11:14
    m=m+1;
    M=rem(m,row_col(2));  M(M==0)=4;
    mm=mm+1;
    %load data
    %load Matt's data and set it in a proper format.
    TrajData=cell(1,1); %cell(10,1); %Here we set the number of realizations (10).
    X_th=X_Th(M);
    ll=0;
    for j=0:(size(TrajData,1)-1)
        ll=ll+1;
        TrajData{ll,1}(:,:)=readmatrix([loadfigdatafolder 'Trajectories\trajectory [x_th, nu, delta]=[' AuxXthStr{M} ', ' AuxDeltaStr{M,1} ', ' AuxDeltaStr{M,2} '] ' num2str(traj(M)-1) '.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<X_th) = Aux(NC<X_th); ND2(NC>=X_th) = Aux(NC>=X_th);
            Kval = k_all{ij}(:,jj);
            nsubrow=3; nsubcol=4;
            ax1(jj)=subplot(nsubrow,nsubcol,jj); hold on;
            plot(Time,X_th*Kval,'b--');
            plot(Time,(1-X_th)*Kval,'r--');
            plot(Time,Kval,'k')
            plot(Time,NC,'b')
            plot(Time,ND1,'r')
            plot(Time,ND2,'g')
        end
    end
    %traj(m)=2;
    ax=ax1(1);
    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]);
    gto_aux2=findobj(gca,'Type','Line',{'LineStyle','--','Color',[0,0,1]});
    gto_aux3=findobj(gca,'Type','Line',{'LineStyle','--','Color',[1,0,0]});
    xlim_aux=200;
    if nu_delta(M,1)<1e2
        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,'LineStyle',':','XData',[XdAta(1),xlim_aux],'YData',CurlyK_aux*[1,1],'LineWidth',2);
        set(gto_aux2,'XData',[XdAta(1),xlim_aux],'YData',X_th*CurlyK_aux*[1,1]);
        set(gto_aux3,'XData',[XdAta(1),xlim_aux],'YData',(1-X_th)*CurlyK_aux*[1,1]);
        set(gca,'YLim',[min(ylim),CurlyK_aux*1.5/1.025]);
    end
    % Aesthetics
    axcp(m).Title.String = ['\{' num2str(round(X_Th(M),2)) ', ' num2str(nu_delta(M,1)) ', ' num2str(nu_delta(M,2)) '\}'];
    xlabel('Time (a.u.)');
    if M==1
        axcp(m).Title.String = ['\{x_{th},\nu,\delta\} = ' axcp(m).Title.String];
        ylabel('Individuals');
    end
    xLim=xlim; xlim([xLim(1) min(xLim(2),xlim_aux)]);
    if M~=row_col(2); yticks([0 kKminus kKplus]); else; yticks([0 kKminus round(CurlyK_aux)]); end
    yticklabels(strsplit(num2str(yticks)));
    if M>1; ylabel(''); end
    set(gca,'Fontsize',axisfontsize);
    set(gca,'YLim',1.025*get(gca,'YLim')); ylim([0 max(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

% Set the legend for some panels
text(axcp(1),0.8,15/100,'x_{th}','Units','Normalized','HorizontalAlignment','center',...
        'VerticalAlignment','bottom','Color',[0,0,0],'FontWeight','Bold','FontSize',panelletterfontsize)
text(axcp(1),0.025,30/100,'$$x=\frac{N_{R}}{N}$$','Units','Normalized','HorizontalAlignment','left',...
        'VerticalAlignment','bottom','Color',[0,0,1],'FontWeight','Bold','FontSize',panelletterfontsize,'interpreter','latex')
text(axcp(5),1.01,100/max(axcp(5).YLim),'K','Units','Normalized','HorizontalAlignment','left',...
        'VerticalAlignment','middle','Color',[0,0,0],'FontWeight','Bold','FontSize',panelletterfontsize)
text(axcp(5),0.5,85/max(axcp(5).YLim),'N','Units','Normalized','HorizontalAlignment','center',...
        'VerticalAlignment','top','Color',[0,0.5,0],'FontWeight','Bold','FontSize',panelletterfontsize)

text(axcp(9),0.55,85/max(axcp(5).YLim),'N_{S}','Units','Normalized','HorizontalAlignment','center',...
        'VerticalAlignment','top','Color',[1,0,0],'FontWeight','Bold','FontSize',panelletterfontsize)
text(axcp(9),0.20,10/max(axcp(5).YLim),'N_{R}','Units','Normalized','HorizontalAlignment','left',...
        'VerticalAlignment','bottom','Color',[0,0,1],'FontWeight','Bold','FontSize',panelletterfontsize)

text(axcp(8),1.01,308/max(axcp(8).YLim),'$\mathcal{K}$','Units','Normalized','HorizontalAlignment','left',...
        'VerticalAlignment','middle','Color',[0,0,0],'FontWeight','Bold','FontSize',panelletterfontsize,'Interpreter','latex')
text(axcp(12),1.01,308/max(axcp(12).YLim),'$\mathcal{K}$','Units','Normalized','HorizontalAlignment','left',...
        'VerticalAlignment','middle','Color',[0,0,0],'FontWeight','Bold','FontSize',panelletterfontsize,'Interpreter','latex')

%% Margins
%set(f,'PaperPositionMode','auto','PaperOrientation','landscape','position',[77.8, 80.2, 1356.8, 380]);
fig_pos = f.PaperPosition;
set(f,'PaperSize',[0.84*fig_pos(3) .97*fig_pos(4)]);
f.PaperPosition(1:2) = f.PaperPosition(1:2) + 0*[.15 .18];

%% Save Figure
print(f,[figfolder figname '_newmargins'],'-dpdf','-painters'); %'-painters' sets de render and saves it as vectorial.
%savefig(f,[figfolder figname '.fig'])
saveas(f,[figfolder figname '.png'],'png')
saveas(f,[figfolder figname '.pdf'],'pdf')
