clearvars
close all
axisfontsize=16;
panelletterfontsize=18;
cd('C:\Users\matlh\OneDrive - University of Leeds\Research\My Projects\In-Leeds\Coop-FixRelThr\Research Data Leeds Repository\')
figfolder= 'C:\Users\matlh\OneDrive - University of Leeds\Research\My Projects\In-Leeds\Coop-FixRelThr\Research Data Leeds Repository\';
%loadfigdatafolder='C:\Users\matlh\OneDrive - University of Leeds\Research\MATLAB\Results\DensDepSys\StaticEnv_ctt_xth\';
runtheory = true;
resultsversion = '_200to25'; %['_' num2str(1)];
savetheoryresults = true;
loadtheory = false;

%% Run Moran theory
if runtheory
    % Set parameteres
    S = 0.1;
    AM = .25;
    Nth = 0:200; %0:400;
    Nmax = max(Nth);
    TheoryData(1).Nc0=0:Nmax; TheoryData(1).Nth=Nth; TheoryData(1).S_AM=cat(1,S,AM); TheoryData(1).N=Nmax*2.^[0 -1 -2 -3]; %Nmax*2.^[0 -1 -2 -3 -4];
    TheoryData(1).FixProbAtThr=nan(numel(Nth),numel(S),numel(TheoryData(1).N)); TheoryData(1).ApproxFixProbAtThr=TheoryData(1).FixProbAtThr; TheoryData(1).ApproxMATAtThr=TheoryData(1).FixProbAtThr; TheoryData(1).MATAtThr=TheoryData(1).FixProbAtThr;
    TheoryData(1).QSCD4FixXthSys=nan(Nmax+1,numel(Nth),numel(S),numel(TheoryData(1).N)); TheoryData(1).ApproxQSCD4FixXthSys=TheoryData(1).QSCD4FixXthSys;
    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);
        n_idx=0;
        for n=TheoryData(1).N
            n_idx=n_idx+1;
            aux=TheoryData(1).Nc0<=n;
            nth_idx=0;
            for nth=TheoryData(1).Nth
                nth_idx=nth_idx+1;
                [TheoryData(1).ApproxFixProbAtThr(nth_idx,s_a_idx,n_idx),TheoryData(1).FixProbAtThr(nth_idx,s_a_idx,n_idx)]=...
                    FixProbAtThr4AsymDensDepSys(n,nth,s,a);
                [TheoryData(1).ApproxMATAtThr(nth_idx,s_a_idx,n_idx),TheoryData(1).MATAtThr(nth_idx,s_a_idx,n_idx),TheoryData(1).SemiExactMATAtThr(nth_idx,s_a_idx,n_idx),TheoryData(1).Approx2MATAtThr(nth_idx,s_a_idx,n_idx)]=...
                    MeanAbsorptionTimeAtThr4AsymDensDepSys(n,nth,s,a);
                [TheoryData(1).ApproxQSCD4FixXthSys(aux,nth_idx,s_a_idx,n_idx),TheoryData(1).QSCD4FixXthSys(aux,nth_idx,s_a_idx,n_idx)]=...
                    QSCD4FixXthSys(TheoryData(1).Nc0(aux),n,nth,s,a);
            end
        end
    end; clearvars nth nth_idx n aux n_idx s a s_a_idx S AM Nth Nmax
    if savetheoryresults
        save([figfolder 'Fig2TheoryData' resultsversion '.mat'],'TheoryData','-v7.3')
    end
elseif loadtheory
    load([figfolder 'Fig2TheoryData' resultsversion '.mat']);
else
    error('Error: data not computed nor loaded.')
end

%(OLD) %% Load Moran theory.
%load([loadfigdatafolder 'AsymDensDepSys_StaticEnv_ctt_xth_Theory_9.mat']);

%% 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 x
    subplot(1,2,1);%subplot(2,2,1)
    hold on;
    LineStyle = {':','--'};
    Col=colormap('summer'); N_Col= size(Col,1); N_n=numel(TheoryData(1).N)+1;
    %Color = {'b','g','r'};
    N_fair = (TheoryData(1).Nth*log(1-a)-log(s*(1-a)/(a-s)))/log(1-s);
    xth_fair = (log(1-s)+(1./TheoryData(1).N).*log(s*(1-a)/(a-s)))/log(1-a);
    Legtxt=cellfun(@num2str,num2cell(TheoryData(1).Nth),'un',0); Legtxt{1}=cat(2,'N_{th}= ',Legtxt{1});
    x_inf_fair = log(1-s)/log(1-a);
    Line(4,numel(TheoryData(1).N),s_a_idx)=plot(x_inf_fair*[1,1],[0 1],'k:','LineWidth',1.5);%Color{nth_idx},'d','LineWidth',2);
    text(x_inf_fair,1,'   $$\frac{\ln{(1-s)}}{\ln{(1-a)}}~$$','Color','k','HorizontalAlignment','right','VerticalAlignment','middle','FontSize',axisfontsize*14/16,'FontWeight','bold', 'Interpreter','latex')
    n_idx=0;
    for n=TheoryData(1).N
        n_idx=n_idx+1;
        Y_approx=TheoryData(1).ApproxFixProbAtThr; Y_approx = Y_approx(:,1,n_idx);
        Y_exact=TheoryData(1).FixProbAtThr; Y_exact = Y_exact(:,1,n_idx);
        nth=TheoryData(1).Nth;
        Line(1,n_idx,s_a_idx)=plot(nth/n,Y_exact,'Color',Col(round(N_Col*(n_idx+1)/(N_n+1)),:),'LineStyle',LineStyle{1},'LineWidth',2);%Color{nth_idx},'LineStyle',LineStyle{2},'LineWidth',2);
        Line(2,n_idx,s_a_idx)=plot(nth/n,Y_approx,'Color',Col(round(N_Col*(n_idx+1)/(N_n+1)),:),'LineStyle',LineStyle{2},'LineWidth',2);%Color{nth_idx},'LineStyle',LineStyle{1},'LineWidth',2);
        Line(3,n_idx,s_a_idx)=scatter(xth_fair(n_idx),.5,50,Col(round(N_Col*(n_idx+1)/(N_n+1)),:),'d','LineWidth',2);%Color{nth_idx},'d','LineWidth',2);
        Line(1,n_idx,s_a_idx).DisplayName=num2str(n); Line(2,n_idx,s_a_idx).DisplayName=num2str(n);
        %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([0 1])%set(gca,'XLim',[min(TheoryData(1).Nth),max(N_fair)+min(N_fair)-min(TheoryData(1).Nth)])
    xticks(0:.25:1)
    xlabel('Relative coop. threshold x_{th}')%xlabel('Total population N')
    ylabel('Resistant fixation prob. \phi')
    ylim(0:1)
    yticks(0:.5:1)
    set(gca,'Fontsize',axisfontsize)
    lg(l)=legend(Line(1,:,1),cellfun(@num2str,num2cell(TheoryData(1).N),'un',0),'box','off','Location','SouthEast');
    title(lg(l),'N')
    lg(1).Visible='on';
    set(Line(3,:,s_a_idx),'Visible','off');
    clearvars Line    
    % Subplot MAT vs x
    subplot(1,2,2);%subplot(2,2,2)
    hold on;
    LineStyle = {':','--'};
    Col=colormap('summer'); N_Col= size(Col,1); N_n=numel(TheoryData(1).N)+1;
    %Color = {'b','g','r'};
    N_fair = (TheoryData(1).Nth*log(1-a)-log(s*(1-a)/(a-s)))/log(1-s);
    xth_fair = (log(1-s)+(1./TheoryData(1).N).*log(s*(1-a)/(a-s)))/log(1-a);
    Legtxt=cellfun(@num2str,num2cell(TheoryData(1).Nth),'un',0); Legtxt{1}=cat(2,'N_{th}= ',Legtxt{1});
    x_inf_fair = log(1-s)/log(1-a);
    Line(4,numel(TheoryData(1).N),s_a_idx)=plot(x_inf_fair*[1,1],[10^0 10^12],'k:','LineWidth',1.5);%Color{nth_idx},'d','LineWidth',2);
    %text(x_inf_fair,10^12,'   $$\frac{\ln{(1-s)}}{\ln{(1-a)}}~$$','Color','k','HorizontalAlignment','right','VerticalAlignment','middle','FontSize',axisfontsize,'FontWeight','bold', 'Interpreter','latex')
    n_idx=0;
    for n=TheoryData(1).N
        n_idx=n_idx+1;
        Y_approx=TheoryData(1).ApproxMATAtThr; Y_approx = Y_approx(:,1,n_idx);
        Y_exact=TheoryData(1).MATAtThr; Y_exact = Y_exact(:,1,n_idx);
        nth=TheoryData(1).Nth;
        Line(1,n_idx,s_a_idx)=plot(nth/n,Y_exact,'Color',Col(round(N_Col*(n_idx+1)/(N_n+1)),:),'LineStyle',LineStyle{1},'LineWidth',2);%Color{nth_idx},'LineStyle',LineStyle{2},'LineWidth',2);
        Line(2,n_idx,s_a_idx)=plot(nth/n,Y_approx,'Color',Col(round(N_Col*(n_idx+1)/(N_n+1)),:),'LineStyle',LineStyle{2},'LineWidth',2);%Color{nth_idx},'LineStyle',LineStyle{1},'LineWidth',2);
        Line(3,n_idx,s_a_idx)=scatter(xth_fair(n_idx),interp1(nth/n,Y_exact,xth_fair(n_idx)),50,Col(round(N_Col*(n_idx+1)/(N_n+1)),:),'d','LineWidth',2);%Color{nth_idx},'d','LineWidth',2);
        Line(1,n_idx,s_a_idx).DisplayName=num2str(n); Line(2,n_idx,s_a_idx).DisplayName=num2str(n);
        %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
    set(gca,'YScale','log')
    xlim([0 1])%set(gca,'XLim',[min(TheoryData(1).Nth),max(N_fair)+min(N_fair)-min(TheoryData(1).Nth)])
    xticks(0:.25:1)
    xlabel('x_{th}')%xlabel('Total population N')
    ylabel('Mean coexistence time (a.u.)')
    ylim(10.^[0 12]); %yticks(0:.5:1);
    set(gca,'Fontsize',axisfontsize)
    lg(l)=legend(Line(1,:,1),cellfun(@num2str,num2cell(TheoryData(1).N),'un',0),'box','off','Location','NorthEast');
    title(lg(l),'N')
    lg(1).Visible='on';
    set(Line(3,:,s_a_idx),'Visible','off');
    clearvars Line    
end
f=gcf;
set(f,'PaperPositionMode','auto','PaperOrientation','landscape','position',[77.8, 80.2, 1356.8, 539.2]);%[77.8, 80.2, 1356.8, 380]);%[13 53 1514.4*2/3 714.4*(2/3)^2]);
Axes=flip(findobj(f,'Type','Axes'));
text(Axes(1),-.05,1.035,'(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.035,'(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';

%% Load and plot Simulaton Data for full model at constant K_0.
Data=readmatrix([figfolder 'SimulationData\fixation probability and time N=[25, 50, 100, 200].csv']); %Data=readmatrix([figfolder 'SimulationData_10xthpoints_1000realisations.csv']);%Data=readmatrix([figfolder 'SimulationData_100realisations.csv']); % load the data
Col=colormap('summer'); N_Col= size(Col,1); N_n=numel(TheoryData(1).N)+1;
n_idx=0;
for n=flip(reshape(unique(Data(:,2)),1,[]))
    n_idx=n_idx+1;
    idx = Data(:,2)==n;
    X = Data(idx,1);
    Y_fixation = Data(idx,3);
    Y_meancoextime = Data(idx,4);
    Color = get(findobj(Axes(1),'Type','Line','DisplayName',num2str(n)),'Color');
    hold on; LIne(1,n_idx)=plot(Axes(1),X,Y_fixation,'-','LineWidth',2,'Color',Col(round(N_Col*(n_idx+1)/(N_n+1)),:)); LIne(1,n_idx).DisplayName=num2str(n);
    plot(Axes(2),X,Y_meancoextime,'-','LineWidth',2,'Color',Col(round(N_Col*(n_idx+1)/(N_n+1)),:));
end
ylim(Axes(2),10.^[0 6])
lg(l)=legend(Axes(1),LIne(1,:),LIne(1,:).DisplayName,'box','off','Location','SouthEast');
title(lg(l),'N')
lg(1).Visible='on';

%% Margins
set(f,'PaperPositionMode','auto','PaperOrientation','landscape','position',[77.8, 80.2, 1356.8, 539.2]); % paper units are centimeters by default.
fig_pos = f.PaperPosition;
set(f,'PaperSize',[0.83*fig_pos(3) 0.99*fig_pos(4)]);
f.PaperPosition(1:2) = f.PaperPosition(1:2) + [.15 .0];

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


%% THEORY
% Functions computing the EXACT and APPROXIMATE fixation probability starting at the threshold in the Asymmetric Density-Dependent System where only when cooperators' fraction is x>xc they generate a ctt public good.
function [y_approx,y_exact] = FixProbAtThr4AsymDensDepSys(N,Nth,s,a)
y_aux = 1-((1-a)/(1-s))^(Nth); % contributions to the cumulative p-distro starting below xc.
y_exact = y_aux/(y_aux + (a-s)/(s*(1-a))*((1-a)/(1-s))^(Nth)*((1/(1-s))^(N-Nth)-1)   ); % contributions to the cumulative p-distro starting above xc.
y_approx = 1/(1+ exp( Nth*log(1-a) - N*log(1-s) + log((a-s)/(s*(1-a))) ) );
end

% Functions computing the APPROXIMATE mean absortion time starting at the threshold in the Asymmetric Density-Dependent System where only when cooperators' fraction is x>xc they generate a ctt public good.
function [y_approx_fixprobpropterm,y_exact,y_semiexact_fixprobpropterm,y_approx_genharmonicseries] = MeanAbsorptionTimeAtThr4AsymDensDepSys(N,Nth,s,a)
[~,FixProb_exact] = FixProbAtThr4AsymDensDepSys(N,Nth,s,a); % compute exact fixation probability from threshold.
y_semiexact_fixprobpropterm = FixProb_exact*SumSumInvTcProdGamaAtThr(N,Nth,s,a);
y_exact = FixProb_exact*SumSumInvTcProdGamaAtThr(N,Nth,s,a)-(1-FixProb_exact)*SumSumInvTcProdGamaNdownNup(1,Nth,N,Nth,s,a); % compute semiexact mean absorption time from threshold.
y_approx_fixprobpropterm = FixProb_exact/s*( (1/(1-s))^(N-Nth) - 1 )*( 1/(N-Nth)*(a*(1-s)/(s*(a-s)) + 0/(N-Nth)*((1-s)^2/s^2-(1-a)*(1-s)/(a-s)^2) ) + 1/Nth*(a*(1-s)/(s*(a-s)) + 0/Nth*((1-a)^2/(a-s)^2-(1-s)/s^2) ) );
c = (1-s)/(1-a);
y_approx_genharmonicseries = FixProb_exact/s*( (1/(1-s))^(N-Nth) - 1 )*( exp((N-Nth+1)*log(c)+ApproxLogGenHarmSeries(1/c,N-Nth+1,N-1,2)) + exp((-Nth)*log(c)+ApproxLogGenHarmSeries(c,1,Nth-1,2)) + exp((N+1-Nth)*log(1-s)+ApproxLogGenHarmSeries(1/(1-s),1,N-Nth,2)) + exp((-Nth)*log(1-s)+ApproxLogGenHarmSeries(1-s,Nth,N-1,2)) );
    function [LogSum] = ApproxLogGenHarmSeries(c,d,u,terms)
        aux = (1:4)<=terms;
        if abs(c)>1
            LogSum = (u+1)*log(c) -log(u*(c-1)) +log( aux(1)*1 + aux(2)*1/(u*(c-1)) + aux(3)*(1+c)/(u*(c-1))^2 + aux(4)*(1+4*c+c^2)/(u*(c-1))^3 );
        elseif abs(c)<1
            LogSum = (d+1)*log(c) -log(d*(1-c)) +log( aux(1)*1/c - aux(2)*1/(d*(1-c)) + aux(3)*(1+c)/(d*(1-c))^2 - aux(4)*(1+4*c+c^2)/(d*(1-c))^3 );
        else
            error('Error: c=1 doesn''t apply to the biophysical model.')
        end
    end
    function [y_approx,y_exact] = FixProbAtThr4AsymDensDepSys(N,Nth,s,a)
        y_aux = 1-((1-a)/(1-s))^(Nth); % contributions to the cumulative p-distro starting below xc.
        y_exact = y_aux/(y_aux + (a-s)/(s*(1-a))*((1-a)/(1-s))^(Nth)*((1/(1-s))^(N-Nth)-1)   ); % contributions to the cumulative p-distro starting above xc.
        y_approx = 1/(1+ exp( Nth*log(1-a) - N*log(1-s) + log((a-s)/(s*(1-a))) ) );
    end
    function T = SumSumInvTcProdGamaAtThr(N,Nth,s,a)
        % computing Sum_{k=Nth}^{N-1}[ Sum_{n=1}^{k}[ Prod_{m=n+1}^{k}[ gama(m,Nth,s) ] /Tc(n) ] ].
        T=0;
        % Terms with n>=Nth in the above expression.
        S1=0;
        for k=Nth:N-1
            S2=0;
            for n=Nth:k
                S2 = S2 + (1-s)^n*( (1-s)/(N-n) + 1/n );
            end
            S1 = S1 + (1/(1-s))^k*S2;
        end
        T = T + 1/(1-s)*S1;
        % Terms with n<Nth in the above expression.
        S1=0;
        for n=1:(Nth-1)
            S1 = S1 + ((1-a)/(1-s))^(-n)*( (1-s)/(N-n) + (1-a)/n );
        end
        S1 = S1 * (1-s)/s*(1-a)^(Nth-1)*((1/(1-s))^N-(1/(1-s))^Nth);
        T = T + 1/(1-s)*S1;
    end
    function T = SumSumInvTcProdGamaNdownNup(Ndown,Nup,N,Nth,s,a)
        % computing Sum_{k=Ndown}^{Nup-1}[ Sum_{n=1}^{k}[ Prod_{m=n+1}^{k}[ gama(m,Nth,s) ] /Tc(n) ] ].
        T=0;
        % Second term in the above expression (notebook) for both below and above Nth.
        S1=0;
        for k=round(((Ndown+Nth)+abs(Ndown-Nth))/2):Nup-1
            S2=0;
            for n=Nth:k
                S2 = S2 + (1-s)^n*( (1-s)/(N-n) + 1/n );
            end
            S1 = S1 + (1/(1-s))^k*S2;
        end
        T = T + 1/(1-s)*S1;
        % Last term in the expression.
        S1=0;
        for n=1:(Nth-1)
            S1 = S1 + ((1-a)/(1-s))^(-n)*( (1-s)/(N-n) + (1-a)/n );
        end
        S1 = S1 * (1-s)/s*(1-a)^(Nth-1)*((1/(1-s))^Nup-(1/(1-s))^(((Ndown+Nth)+abs(Ndown-Nth))/2));
        T = T + 1/(1-s)*S1;
        % First term in the expression that only appears below Nth.
        S1=0;
        for k=Ndown:Nth-1
            S2=0;
            for n=1:k
                S2 = S2 + ((1-a)/(1-s))^(-n)*( (1-s)/(N-n) + (1-a)/n );
            end
            S1 = S1 + ((1-a)/(1-s))^k*S2;
        end
        T = T + 1/(1-s)*S1;
    end
end

% Functions computing the EXACT and APPROXIMATE Quasi-stationary Composition Distribution the fixed threshold case, starting at the threshold in the Asymmetric Density-Dependent System where only when cooperators' fraction is x>xc they generate a ctt public good.
function [y_approx,y_exact] = QSCD4FixXthSys(Nc,N,nth,s,a)
Nc((Nc==0)|(Nc==N))=nan;
aux1 = ((1-s)/(1-a)).^(((nth-1+Nc)-abs(nth-1-Nc))/2);
beta = ( (nth-1)*(N-nth+1)*(N-s*nth) )/( ( (1-a)*N+(a-s)*(nth-1) )*(nth)*(N-nth) );
beta_approx = (1-s*nth/N)/((1-a)+(a-s)*nth/N);%(1-s*nth/N)/(1-a);
aux2 = beta.^(Nc>=nth);
aux2_approx = beta_approx.^(Nc>=nth);%beta.^(Nc>=nth);
aux3= (1-s).^(((-nth+1+Nc)+abs(-nth+1+Nc))/2);
y_approx = aux1.*aux2_approx.*aux3/( (s+(a-s)*beta_approx)*(1-s)/(s*(a-s))*((1-s)/(1-a))^(nth-1) );%aux1.*aux2_approx.*aux3/( (s+(a-s)*beta)*(1-s)/(s*(a-s))*((1-s)/(1-a))^(nth-1) );
y_exact = aux1.*aux2.*aux3/( (1-s)/(a-s)*(((1-s)/(1-a))^(nth-1)-1) + (1-s)/s*((1-s)/(1-a))^(nth-1)*beta*(1-(1-s)^(N-nth)));
end