clearvars
close all
axisfontsize=16;
panelletterfontsize=18;
figfolder= 'C:\Users\matlh\OneDrive - University of Leeds\My Projects\In-Leeds\Cooperator-Treshold system\Research Data Leeds Repository\';
loadfigdatafolder='C:\Users\matlh\OneDrive - University of Leeds\My Projects\In-Leeds\Cooperator-Treshold system\Research Data Leeds Repository\SimulationData\';
nrealizations=1000;

%% Load Moran theory.
load([loadfigdatafolder 'FluctNcttK\AsymDensDepSys_Moran_Theory.mat']);
TheoryData(1).Nth = TheoryData(1).Nth(1:end-2); %Discard N_th=[80 100] because simulations take too long.

%% Plot FixProb&MAT vs N.
s_a_idx=0;
for s=TheoryData(1).S_AM(1,:)
    s_a_idx=s_a_idx+1;
    a=TheoryData(1).S_AM(2,s_a_idx);
    l=1;
    f(l)=figure; set(gcf,'position',[13 100 1514.4 525],'Name',['AsymDensDepSys_Moran_Theory_FluctNcttK_Simulations_s_' num2str(TheoryData(1).S_AM(1,s_a_idx)) '_a_' num2str(TheoryData(1).S_AM(2,s_a_idx)) '_FixProbAndAbsTimeVsN'])
    % Subplot FixProb vs N
    subplot(1,2,1)
    hold on;
    LineStyle = {'-',':'};
    Col=colormap('summer'); N_Col= size(Col,1); N_nth=numel(TheoryData(1).Nth);
    %Color = {'b','g','r'};
    N_fair = (TheoryData(1).Nth*log(1-a)-log(s*(1-a)/(a-s)))/log(1-s);
    Legtxt=cellfun(@num2str,num2cell(TheoryData(1).Nth),'un',0); Legtxt{1}=cat(2,'N_{th}= ',Legtxt{1});
    x_inf_fair = log(1-a)/log(1-s);
    Line(4,numel(TheoryData(1).Nth),s_a_idx)=plot(x_inf_fair*[1,1],[0 1],'k','LineWidth',2);%Color{nth_idx},'d','LineWidth',2);
    text(x_inf_fair,1,'   $$~~~\frac{\ln{(1-a)}}{\ln{(1-s)}}$$','Color','k','HorizontalAlignment','left','VerticalAlignment','middle','FontSize',axisfontsize,'FontWeight','bold', 'Interpreter','latex')
    nth_idx=0;
    for nth=TheoryData(1).Nth
        nth_idx=nth_idx+1;
        Y_approx=permute(TheoryData(1).ApproxFixProbAtThr,[3 2 1]); Y_approx = Y_approx(:,1,nth_idx);
        Y_exact=permute(TheoryData(1).FixProbAtThr,[3 2 1]); Y_exact = Y_exact(:,1,nth_idx);
        X=permute(TheoryData(1).N,[2 1]);
        Line(1,nth_idx,s_a_idx)=plot(X/nth,Y_approx,'Color',Col(round(N_Col*nth_idx/(N_nth+1)),:),'LineStyle',LineStyle{1},'LineWidth',2);%Color{nth_idx},'LineStyle',LineStyle{1},'LineWidth',2);
        Line(2,nth_idx,s_a_idx)=plot(X/nth,Y_exact,'Color',Col(round(N_Col*nth_idx/(N_nth+1)),:),'LineStyle',LineStyle{2},'LineWidth',2);%Color{nth_idx},'LineStyle',LineStyle{2},'LineWidth',2);
        Line(3,nth_idx,s_a_idx)=scatter(N_fair(nth_idx)/nth,.5,50,Col(round(N_Col*nth_idx/(N_nth+1)),:),'d','LineWidth',2);%Color{nth_idx},'d','LineWidth',2);
        Line(1,nth_idx,s_a_idx).DisplayName=num2str(nth); Line(2,nth_idx,s_a_idx).DisplayName=num2str(nth);
        %text((N_fair(nth_idx)-6)/nth,.75,Legtxt{nth_idx},'Color',Col(round(N_Col*nth_idx/(N_nth+1)),:),'HorizontalAlignment','right','VerticalAlignment','top','FontSize',axisfontsize,'FontWeight','bold')
    end
    xlim([1 4])%set(gca,'XLim',[min(TheoryData(1).Nth),max(N_fair)+min(N_fair)-min(TheoryData(1).Nth)])
    xticks(1:4)
    xlabel('K_{0}/N_{th}')%xlabel('Total population N')
    ylabel('\phi')
    ylim(0:1)
    yticks(0:.5:1)
    set(gca,'Fontsize',axisfontsize)
    %lg(l)=legend(Line(2,:,1),cellfun(@num2str,num2cell(TheoryData(1).Nth),'un',0),'box','off','Location','SouthWest');
    %title(lg(l),'N_{th}')
    %lg(1).Visible='on';
    clearvars Line
    % Subplot MAT vs N
    subplot(1,2,2)
    hold on;
    LineStyle = {'-',':','-.'};
    Col=colormap('summer'); N_Col= size(Col,1); N_nth=numel(TheoryData(1).Nth);
    %Color = {'b','g','r'};
    N_fair = (TheoryData(1).Nth*log(1-a)-log(s*(1-a)/(a-s)))/log(1-s);
    N_maxt = N_fair+log(TheoryData(1).Nth)/log(1/(1-s));
    x_inf_fair = log(1-a)/log(1-s);
    Line(3,nth_idx,s_a_idx)=plot(x_inf_fair*[1,1],[1 10^(20)],'k','LineWidth',2);%Color{nth_idx},'d','LineWidth',2);
    %text(x_inf_fair,10^(20),'$$\frac{K_{0}^*}{N_{th}}$$~~~','Color','k','HorizontalAlignment','right','VerticalAlignment','middle','FontSize',axisfontsize,'FontWeight','bold', 'Interpreter','latex')
    nth_idx=0;
    for nth=TheoryData(1).Nth
        nth_idx=nth_idx+1;
        Y_exact=permute(TheoryData(1).MATAtThr,[3 2 1]); Y_exact = Y_exact(:,1,nth_idx);
        %%Y_semiexact_proptofixprob=permute(TheoryData(1).SemiExactMATAtThr,[3 2 1]); Y_semiexact_proptofixprob = Y_semiexact_proptofixprob(:,1,nth_idx);
        Y_approx=permute(TheoryData(1).ApproxMATAtThr,[3 2 1]); Y_approx = Y_approx(:,1,nth_idx);
        Y_approx2=permute(TheoryData(1).Approx2MATAtThr,[3 2 1]); Y_approx2 = Y_approx2(:,1,nth_idx);
        X=permute(TheoryData(1).N,[2 1]);
        Line(1,nth_idx,s_a_idx)=plot(X/nth,Y_exact,'Color',Col(round(N_Col*nth_idx/(N_nth+1)),:),'LineStyle',LineStyle{1},'LineWidth',2);%Color{nth_idx},'LineStyle',LineStyle{1},'LineWidth',2);
        Line(2,nth_idx,s_a_idx)=scatter(N_fair(nth_idx)/nth,interp1(X,Y_exact,N_fair(nth_idx)),50,Col(round(N_Col*nth_idx/(N_nth+1)),:),'d','LineWidth',2);%Color{nth_idx},'d','LineWidth',2);
        %Line(2,nth_idx,s_a_idx)=scatter(N_maxt(nth_idx),interp1(X,Y_exact,N_maxt(nth_idx)),50,Col(round(N_Col*nth_idx/(N_nth+1)),:),'^','LineWidth',2);%Color{nth_idx},'^','LineWidth',2);
        %%plot(X,Y_semiexact_proptofixprob,Color{nth_idx},'LineStyle',LineStyle{nth_idx},'LineWidth',2)
        %Line(4,nth_idx,s_a_idx)=plot(X/nth,Y_approx,'color',Col(round(N_Col*nth_idx/(N_nth+1)),:),'LineStyle',LineStyle{2},'LineWidth',1);
        %Line(5,nth_idx,s_a_idx)=plot(X,Y_approx2,Color{nth_idx},'LineStyle',LineStyle{3},'LineWidth',1);
    end
    set(gca,'YScale','log')
    xlim([1 4]);%xlim([min(TheoryData(1).Nth),max(N_maxt)+1.5*(min(N_maxt)-.95*max(TheoryData(1).Nth))])
    xticks(1:4)
    %ylim([10^(-1) 10^(7)])
    xlabel('K_{0}/N_{th}')%xlabel('Total population N')
    ylabel('MCT (a.u.)')%ylabel('Mean cabsorption time')
    set(gca,'Fontsize',axisfontsize)
    lg(l)=legend(Line(1,:,1),cellfun(@num2str,num2cell(TheoryData(1).Nth),'un',0),'box','off','Location','NorthWest');
    title(lg(l),'N_{th}')
    lg(l).Visible='off';
    set(gca,'YScale','log')
    ax(1)=gca;
    ax(2)=axes('Outerposition',[.57 0.62 .15 .3]);%ax(2)=axes('Outerposition',[.65 0.175 .275 .3]);
    box on
    copyobj(ax(1).Children,ax(2));
    set(gca,'YScale','lin')
    xlim([2 4]);%xlim([min(TheoryData(1).Nth),max(TheoryData(1).N)])%xlim([min(TheoryData(1).Nth),max(N_maxt)+1.5*(min(N_maxt)-.95*max(TheoryData(1).Nth))])
    ylim([0 5*10^(19)])
end
f=gcf;
set(f,'PaperPositionMode','auto','PaperOrientation','landscape','position',[13 53 .875*1514.4 714.4*7/10]);
Axes=flip(findobj(f,'Type','Axes'));
text(Axes(1),-.05,1.05,'(a)','Units','Normalized','HorizontalAlignment','right',...
    'VerticalAlignment','bottom','Color',[0,0,0],'FontWeight','Bold','FontSize',panelletterfontsize)
set(Axes(1),'Fontsize',axisfontsize)
text(Axes(2),-.05,1.05,'(b)','Units','Normalized','HorizontalAlignment','right',...
    'VerticalAlignment','bottom','Color',[0,0,0],'FontWeight','Bold','FontSize',panelletterfontsize)
set(Axes(2),'Fontsize',axisfontsize)
Axes(2).Legend.Visible = 'off';
arrayfun(@(x)xlabel(x,'K_{0}/N_{th}'),Axes(1:2),'un',0);
Axes(1).YTick=0:.2:1; Axes(1).XTick=1:4;
Axes(2).YLabel.String='MCT (a.u.)'; Axes(2).YLim=[1 10^(9)]; Axes(2).XTick=1:4;
delete(Axes(3)); % Axes(3).XLim=[35 1500]; Axes(3).XScale='log'; %Axes(3).XTick=1.5*[10,100,1000]; Axes(3).XTickLabel=cellfun(@num2str,num2cell(Axes(3).XTick),'un',0);
aux=findobj(Axes(1:3),'Type','Scatter'); delete(aux);
aux=findobj(Axes(1),'Type','Line','LineStyle','-',{'-not','Color','k'}); delete(aux);
aux=findobj(Axes(2),'Type','Line','LineStyle','-',{'-not','Color','k'}); set(aux(:),'Linestyle',':');

%% Load and Preprocess Simulaton Data for full model at constant K_0.
Data=nan(364-1,4,1000);
tic
for i=1:nrealizations
    Data(:,:,i)=readmatrix([loadfigdatafolder 'FluctNcttK\fixation probability and time ' num2str(i) '.csv']);
end
Data=mean(Data,3);
disp(['Time to load and preprocess simulation data' num2str(toc)])

%% Plot the data
N_th=reshape(unique(Data(:,1)),1,[]);
nth_idx=0;
for nth=N_th
    nth_idx=nth_idx+1;
    idx = abs(Data(:,1)-nth)<1e-6;
    X = Data(idx,2)./Data(idx,1);
    Y_fixation = Data(idx,3);
    Y_meancoextime = Data(idx,4);
    Color = get(findobj(Axes(1),'Type','Line','DisplayName',num2str(nth)),'Color');
    hold on; LIne(1,nth_idx)=plot(Axes(1),X,Y_fixation,'-','LineWidth',2,'Color',Color); LIne(1,nth_idx).DisplayName=num2str(nth);
    plot(Axes(2),X,Y_meancoextime,'-','LineWidth',2,'Color',Color);
end
ylim(Axes(2),10.^[0 6])
lg(l)=legend(Axes(1),LIne(1,:),LIne(1,:).DisplayName,'box','off','Location','SouthWest');
title(lg(l),'N_{th}')
lg(1).Visible='on';
    
%% Save Figure
savefig(f,[figfolder 'FinalFigS1.fig'])
saveas(f,[figfolder 'FinalFigS1.png'],'png')
saveas(f,[figfolder 'FinalFigS1.pdf'],'pdf')