clearvars
%close all
colormapfontsize=14;
axisfontsize=16;
panelletterfontsize=18;
runtheory = false;%true;%false;
savetheoryresults = false;%true;%false;
loadtheory = true;%true;%false;
nu_rescaling = false;%true;%false; % rescaling with nu->nu/(a-s) for xth<ln(1-s)/ln(1-a) and nu->nu/s otherwise.
K.values = [100; 1000];%[25; 200];%[25; 400];%[100; 1000];
resultsversion = ['_K_' num2str(K.values(1)) '_' num2str(K.values(2))]; %['_' num2str(1)];
nu_delta=[1e-4,-.5; 1e-1,-.3; 2e0,.5; 1e2,.7].*[1,0];
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)'};
nrows = 2;%3;
ncols = 2;
s=.1; a=.25; kKminus=K.values(1); kKplus=K.values(2);
figidx = {12,13,14, 15};
X_TH=(3+[0 1 2 3]*8)/30;
N_TH=repmat(round(X_TH,2),1,nrows); 
aux=cellfun(@(x)cat(2,num2str(x),'_'),num2cell(N_TH),'un',0); aux = [aux{:}]; aux = aux(1:end-1); %THIS DOESN'T MATTER!
Delta = [-1:1]/2;
figname=['FinalFig5'];
loadfigdatafolder='C:\Users\matlh\OneDrive - University of Leeds\Research\My Projects\In-Leeds\Coop-FixRelThr\Research Data Leeds Repository\SimulationData\';
% Initialize figure
f = figure('position',[20 20 1050 2*426.4],'Name',figname); hold on; %figure('position',[13 53 1325.6 nrows/1*595.2],'Name',figname); hold on;
set(f,'PaperPositionMode','auto');
set(f,'PaperOrientation','landscape');
MatrixPos=cat(3,[1,0,.425,0;0,1,0,.65;0,0,.5,0;0,0,0,.25]' , [1,0,.35,0;0,1,0,.25;0,0,.5,0;0,0,0,.25]' , [1,0,1-.5+.15,0;0,1,0,.075;0,0,.5,0;0,0,0,.25]' , [1,0,.45,0;0,1,0,.5;0,0,.5,0;0,0,0,.25]' );
YlIm = [ [0,.01] ; [0,1] ; [.99,1] ; [.99,1] ];
for k=1:nrows*ncols
    ax0(k) = subplot(nrows,ncols,k,'parent',f); % create subplot axes.
    ax0(k).Position = ax0(k).Position + [1-2*rem(k,ncols),(1-2*fix(k/(ncols+1)))/2,0,0]/60; % position them.
    ax0(nrows*ncols+k) = copyobj(ax0(k),f); % copy them to create insets.
    ax0(nrows*ncols+k).Position = ax0(nrows*ncols+k).Position*MatrixPos(:,:,k); % position the insets.
end


titl = {['x_{th} = ' num2str(N_TH(1))], [num2str(N_TH(2))], [num2str(N_TH(3))], [num2str(N_TH(4))]};

%% Compute the theoretical diagrams
if runtheory
    % Set parameteres
    s = 0.1;
    a = .25;
    K.nu = 10.^(-4:(2-(-4))/(101-1):2); K.nu = repmat(K.nu,1,numel(K.nu));
    K.delta = -1:(1-(-1))/(101-1):1; K.delta = reshape(K.delta.*ones(length(K.delta),1),1,[]);
    K2=K;
    %K2.nu = nu_delta(:,1)'; K2.delta = nu_delta(:,2)'; % if we only want "size(nu_delta,1)" dynamic environments.
    Nth = 0:K.values(2); %this sets all the maximum range.
    xthpoints = X_TH;
    TheoryData(1).xth=xthpoints; TheoryData(1).nu_delta=[K.nu;K.delta]'; TheoryData(1).S_AM=cat(1,s,a); TheoryData(1).K=K;
    TheoryData(1).ResAndFixProbAtThr=nan(numel(xthpoints),numel(K2.nu)); TheoryData(1).SensAndFixProbAtThr=TheoryData(1).ResAndFixProbAtThr;
    TheoryData(1).MeanFixRate=nan(numel(xthpoints),numel(K2.nu)); TheoryData(1).MeanExtRate=TheoryData(1).MeanFixRate;
    TheoryData(1).AnnealedSensAndFixNonNormProbArThr=TheoryData(1).MeanFixRate; TheoryData(1).AnnealedResAndFixNonNormProbArThr=TheoryData(1).MeanFixRate;
    TheoryData(1).MeanN=nan(1,numel(K2.nu));
    Xth = []; Aux1 = []; Aux2 = []; Aux3 = []; Aux4 = []; Aux5 = []; Aux6 = [];
    tic
    parfor IDXxth=1:numel(TheoryData(1).xth)
        xth=X_TH(IDXxth);
        Xth = [Xth;xth];
        %aux1=nan(1,numel(K2.nu));aux2=nan(1,numel(K2.nu));
        K3=K2; % prepare the rescaling.
        if nu_rescaling
            if xth<log(1-s)/log(1-a) % for higher sensitive fixation probability.
                K3.nu = K3.nu/(a-s); % rescale nu as in Mauro's PRL2017 using (a-s).
            else
                K3.nu = K3.nu/s; % rescale nu as in Mauro's PRL2017 using s.
            end
        end
        %[aux1,aux2] = NPDMPmarginalization(@FixAndCoexDiagram,xth,s,a,K3,true);%K,true);
        [~,sensfixprob_k]=FixProbAtThr4AsymDensDepSys(kKminus,xth,s,a); resfixprob_k=sensfixprob_k; sensfixprob_k=1-resfixprob_k; [~,mct_k,~,~] = MeanAbsorptionTimeAtThr4AsymDensDepSys(kKminus,xth,s,a); % computations for Kminus static environment.
        [~,sensfixprob_K]=FixProbAtThr4AsymDensDepSys(kKplus,xth,s,a); resfixprob_K=sensfixprob_K; sensfixprob_K=1-resfixprob_K; [~,mct_K,~,~] = MeanAbsorptionTimeAtThr4AsymDensDepSys(kKplus,xth,s,a); % computations for Kplus static environment.
        aux1=(1-K3.delta)/2*sensfixprob_k.*(1-exp(-((1-K3.delta)*kKminus+(1+K3.delta)*kKplus)/mct_k))+(1+K3.delta)/2*sensfixprob_K.*(1-exp(-((1-K3.delta)*kKminus+(1+K3.delta)*kKplus)/mct_K));
        aux2=(1-K3.delta)/2*resfixprob_k.*(1-exp(-((1-K3.delta)*kKminus+(1+K3.delta)*kKplus)/mct_k))+(1+K3.delta)/2*resfixprob_K.*(1-exp(-((1-K3.delta)*kKminus+(1+K3.delta)*kKplus)/mct_K)); 
        Aux1 = [Aux1;aux1];
        Aux2 = [Aux2;aux2];
        [aux3,aux4] = NPDMPmarginalization(@ExtFixRates4FixXthSys,xth,s,a,K3,true);%K,true);
        Aux3 = [Aux3;aux3];
        Aux4 = [Aux4;aux4];
        [aux5,aux6] = NPDMPmarginalization(@AnnealedFixAndCoexDiagram,xth,s,a,K3,true);%K,true);
        Aux5 = [Aux5;aux5];
        Aux6 = [Aux6;aux6];
    end
    [~,iDxx]=sort(Xth);
    TheoryData(1).SensAndFixProbAtThr(iDxx,:) = Aux1; TheoryData(1).ResAndFixProbAtThr(iDxx,:) = Aux2;
    TheoryData(1).MeanFixRate(iDxx,:) = Aux4; TheoryData(1).MeanExtRate(iDxx,:) = Aux3;
    TheoryData(1).AnnealedSensAndFixNonNormProbArThr = Aux5; TheoryData(1).AnnealedResAndFixNonNormProbArThr(iDxx,:) = Aux6;
    clearvars xth Xth sensfixprob_k aUx2_k aUx1_K aUx2_K aux1 Aux1 aux2 Aux2 aux3 Aux3 aux4 Aux4 aux5 Aux5 aux6 Aux6
    toc
    tic
    [TheoryData(1).MeanN,~]=NPDMPmarginalization(@(w,x,y,z)deal(w,0),[],s,a,K2,true); % compute the mean total population for each environment.
    if nu_rescaling
        K3=K2; K3.nu = K3.nu/(a-s); [TheoryData(1).MeanNlowXth,~]=NPDMPmarginalization(@(w,x,y,z)deal(w,0),[],s,a,K3,true); % compute the mean total population for each environment.
        K3=K2; K3.nu = K3.nu/s; [TheoryData(1).MeanNhighXth,~]=NPDMPmarginalization(@(w,x,y,z)deal(w,0),[],s,a,K3,true); % compute the mean total population for each environment.
    end
    %TheoryData(1).MeanN(K2.delta==-1)=K2.values(1); TheoryData(1).MeanN(K2.delta==1)=K2.values(2);
    toc
    clearvars nth nth_idx aux s_a_idx S AM Nth Nmax K2
    if savetheoryresults
        save([figfolder 'Fig4and5TheoryData' resultsversion '.mat'],'TheoryData','-v7.3')
    end
elseif loadtheory
    load([figfolder 'Fig4and5TheoryData' resultsversion '.mat']);
else
    error('Error: data not computed nor loaded.')
end

%% Plot the conditional R fixation probability when there is no coexistence.
m=0;
Clrs(1:3,:) = [ 253 127 57; 0 0 0; 154 205 50]/255;
Clrs2(1:4,:) = [ 1 0 0; .75 0 .75; .25 .25 1; 0 0 .5];
Mkrs = {'v','o','^'}; %Mkrs = {'.','x','*'};
LinSytl = {':','-','--'};
set(ax0(1:2*nrows*ncols),'XScale','log');
nu_npoints=101; delta_npoints=nu_npoints;
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);
kk=0;
for idx_xth=1:numel(X_TH)
    kk=kk+1;
    m=m+1;
    % Theory
    Color=[TheoryData(1).SensAndFixProbAtThr(idx_xth,:)',zeros(size(TheoryData(1).ResAndFixProbAtThr(idx_xth,:)')),TheoryData(1).ResAndFixProbAtThr(idx_xth,:)'];
    MeanFixRate=TheoryData(1).MeanFixRate(idx_xth,:)'; MeanExtRate=TheoryData(1).MeanExtRate(idx_xth,:)'; MeanN=TheoryData(1).MeanN(1,:)';
    AnnealedSensAndFixNonNormProbArThr = TheoryData(1).AnnealedSensAndFixNonNormProbArThr(idx_xth,:)'; AnnealedResAndFixNonNormProbArThr=TheoryData(1).AnnealedResAndFixNonNormProbArThr(idx_xth,:)';
    if nu_rescaling
        if N_TH(idx_xth)<log(1-s)/log(1-a) % for higher sensitive fixation probability.
            MeanN=TheoryData(1).MeanNlowXth(1,:)'; % rescale nu as in Mauro's PRL2017 using nu/(a-s).
        else
            MeanN=TheoryData(1).MeanNhighXth(1,:)'; % rescale nu as in Mauro's PRL2017 using nu/s.
        end
    end
    Color2=(1-exp(-(2*MeanN).*(MeanFixRate+MeanExtRate))).*[MeanExtRate./(MeanFixRate+MeanExtRate),zeros(size(MeanFixRate)),MeanFixRate./(MeanFixRate+MeanExtRate)];
    %Color2=[AnnealedSensAndFixNonNormProbArThr./(MeanFixRate+MeanExtRate),zeros(size(MeanFixRate)),AnnealedResAndFixNonNormProbArThr./(MeanFixRate+MeanExtRate)];
    %Color3=exp(-((1-Data.delta)*kKminus+(1+Data.delta)*kKplus).*Data.nu.*(1-abs(Data.delta))).*Color+(1-exp(-((1-Data.delta)*kKminus+(1+Data.delta)*kKplus).*Data.nu.*(1-abs(Data.delta)))).*Color2;
    Color3=exp(-(2*MeanN).*Data.nu.*(1-abs(Data.delta))).*Color+(1-exp(-(2*MeanN).*Data.nu.*(1-abs(Data.delta)))).*Color2;
    CoexProb1=1-sum(Color(:,[1 3]),2); CoexProb2=1-sum(Color2(:,[1 3]),2);
    %ProbCoex=exp(-((1-Data.delta)*kKminus+(1+Data.delta)*kKplus).*Data.nu.*(1-abs(Data.delta))).*CoexProb1+(1-exp(-((1-Data.delta)*kKminus+(1+Data.delta)*kKplus).*Data.nu.*(1-abs(Data.delta)))).*CoexProb2;
    ProbCoex=exp(-(2*MeanN).*Data.nu.*(1-abs(Data.delta))).*CoexProb1+(1-exp(-(2*MeanN).*Data.nu.*(1-abs(Data.delta)))).*CoexProb2;
    k=0;
    for delta=Delta
        k=k+1;
        delta_idx=Data.delta==delta;
        X=Data.nu(delta_idx); Y=ProbCoex(delta_idx); %log(ProbCoex(delta_idx)./(1-ProbCoex(delta_idx)));
        Y2=Color3(delta_idx,3)./(Color3(delta_idx,1)+Color3(delta_idx,3)); %Y=log(Color3(delta_idx,3)./Color3(delta_idx,1)); %Y=Color3(delta_idx,3)./(Color3(delta_idx,1)+Color3(delta_idx,3));
        axes(ax0(m)); hold on; plot(X,Y,LinSytl{k},'Color',Clrs2(kk,:),'LineWidth',2); %,Clrs(k,:),'LineWidth',2); % plot in panel b.
        axes(ax0(m+nrows*ncols)); hold on; plot(X,Y2,LinSytl{k},'Color',Clrs2(kk,:),'LineWidth',2); %Clrs(k,:),'LineWidth',2); % plot in panel a.
    end
    % Simulations
    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\MATLAB\heatmap K_0=550, gamma=0.82, xth_values=[0.1, 0.36666666666666664, 0.6333333333333333, 0.9] ' num2str(j) '.mat'],'data'));
        %Aux = squeeze(load([loadfigdatafolder 'FixAndCoexProbData\heatmap K_0=550, gamma=0.82, xth_values=[0.1, 0.37, 0.63, 0.9] ' num2str(j) '.mat'],'data'));
        dAta(:,:,:,j) = squeeze(Aux.data(1,idx_xth,:,:,:)); %first dimension is (s,a) parameters, second and onwards as in "data" variable.
        clearvars Aux
    end
    Data.N_C=cell(length(Data.nu),1); Data.N_D=Data.N_C; Data.t=Data.N_C; %Set data in a proper format. Matt's simulations always ran up to t=2<K(nu=0,delta)> even with prior fixation.
    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
    k=0;
    for delta=Delta
        k=k+1;
        delta_idx=Data.delta==delta;
        X=Data.nu(delta_idx);
        N_C=Data.N_C(delta_idx); N_D=Data.N_D(delta_idx); t=Data.t(delta_idx);
        maxt=2*(kKminus*(1-Data.delta(delta_idx))/2+kKplus*(1+Data.delta(delta_idx))/2); %Process data.
        N=cellfun(@(x,y)x+y,N_C,N_D,'un',0);
        NoCoex=cellfun(@(x,y,z)x<min(2*mean(y),max(x)),t,N,num2cell(maxt),'un',0); %cellfun(@(x,y,z)x<min(2*mean(y),z),Data.t,Data.N,num2cell(maxt),'un',0);
        ProbCoex=cellfun(@(x)1-sum(x)/numel(x),NoCoex,'un',1);
        ProbResFixNoCoex=cellfun(@(x,y,z)mean((x(z)+0*1)./(x(z)+y(z)+0*2)),N_C,N_D,NoCoex,'un',1); %cellfun(@(x,y,z)log(mean((x(z)+0*1)./(x(z)+y(z)+0*2))./(1-mean((x(z)+0*1)./(x(z)+y(z)+0*2)))),N_C,N_D,NoCoex,'un',1); %cellfun(@(x,y,z)mean((x(z)+0*1)./(x(z)+y(z)+0*2)),N_C,N_D,NoCoex,'un',1);
        Y=ProbCoex; %log(ProbCoex./(1-ProbCoex));
        axes(ax0(m)); hold on; scatter(X,Y,25,Clrs2(kk,:),Mkrs{k},'LineWidth',1,'MarkerEdgeAlpha',.33); % plot in panel b.
        axes(ax0(m+nrows*ncols)); hold on; scatter(X,ProbResFixNoCoex,25,Clrs2(kk,:),Mkrs{k},'LineWidth',1,'MarkerEdgeAlpha',.33); % plot in panel a.
    end
    axes(ax0(m)); xlabel('\nu'); if rem(m,ncols)==1; ylabel('$$P\left(t>2\langle N\rangle\right)$$','Interpreter','latex'); end
    xlim(10.^[-4,2].*10.^[-.5,.5]); xticks(10.^(-4:2:2)); ylim(0:1);
    title(''); ax0(m).Title.String = titl{m}; ax0(m).Title.Visible = 'on';
    set(gca,'Fontsize',axisfontsize);
    text(ax0(m),-.05,1.05,panelletter{m},'Units','Normalized','HorizontalAlignment','right',...
        'VerticalAlignment','bottom','Color',[0,0,0],'FontWeight','Bold','FontSize',panelletterfontsize)
    chld=findobj(ax0(m),'Type','Line');
    if m==1; LgdStr=flip(cellfun(@num2str,num2cell(Delta),'un',0)); LgdStr{1}=cat(2,'\delta = ',LgdStr{1}); l=legend(ax0(m),chld,LgdStr,'Location','northeast'); l.Box='off'; l.Position=l.Position+[-.1 -.2 0 0]; end %l.Title.String = '\delta';
    if m==nrows*(ncols-1)+1; xlabel('Env. switching freq. \nu (a.u.)'); end
    axes(ax0(m+nrows*ncols)); xlabel('');
    xlim(10.^[-4,2].*10.^[-.5,.5]); xticks(10.^(-4:2:2));
    ylim(YlIm(m,:)); yticks(YlIm(m,:));
    if m==1; ylabel('$$P\left(R|t<2\langle N\rangle\right)$$','Interpreter','latex'); xlabel('\nu'); else; xticklabels({}); end
    set(ax0(m+nrows*ncols),'Fontsize',colormapfontsize);
end

%% Margins
%set(f,'PaperPositionMode','auto','PaperOrientation','landscape','position',[77.8, 80.2, 1356.8, 380]);
fig_pos = f.PaperPosition;
set(f,'PaperSize',[0.88*fig_pos(3) 0.97*fig_pos(4)]); %the latter numbers for the 2-row version.
f.PaperPosition(1:2) = f.PaperPosition(1:2) + [.35 0]; %the latter numbers for the 2-row version.

%% Save Figure
str='';
if nu_rescaling
    str=cat(2,str,'_RescaledNu');
end
print(f,[figfolder figname str],'-dpdf','-painters'); %'-painters' sets de render and saves it as vectorial.
savefig(f,[figfolder figname str '.fig'])
saveas(f,[figfolder figname str '.png'],'png')
%saveas(f,[figfolder figname '.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,xth,s,a)
y_aux = 1-((1-a)/(1-s))^(N*xth); % contributions to the cumulative p-distro starting below xc.
y_exact = y_aux/(y_aux + (a-s)/(s*(1-a))*((1-a)/(1-s))^(N*xth)*((1/(1-s))^(N-N*xth)-1)   ); % contributions to the cumulative p-distro starting above xc.
y_approx = 1/(1+ exp( N*xth*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,xth,s,a)
[~,FixProb_exact] = FixProbAtThr4AsymDensDepSys(N,xth,s,a); % compute exact fixation probability from threshold.
y_semiexact_fixprobpropterm = FixProb_exact*SumSumInvTcProdGamaAtThr(N,N*xth,s,a);
y_exact = FixProb_exact*SumSumInvTcProdGamaAtThr(N,N*xth,s,a)-(1-FixProb_exact)*SumSumInvTcProdGamaNdownNup(1,N*xth,N,N*xth,s,a); % compute semiexact mean absorption time from threshold.
y_approx_fixprobpropterm = FixProb_exact/s*( (1/(1-s))^(N-N*xth) - 1 )*( 1/(N-N*xth)*(a*(1-s)/(s*(a-s)) + 0/(N-N*xth)*((1-s)^2/s^2-(1-a)*(1-s)/(a-s)^2) ) + 1/(N*xth)*(a*(1-s)/(s*(a-s)) + 0/(N*xth)*((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-N*xth) - 1 )*( exp((N-N*xth+1)*log(c)+ApproxLogGenHarmSeries(1/c,N-N*xth+1,N-1,2)) + exp((-N*xth)*log(c)+ApproxLogGenHarmSeries(c,1,N*xth-1,2)) + exp((N+1-N*xth)*log(1-s)+ApproxLogGenHarmSeries(1/(1-s),1,N-N*xth,2)) + exp((-N*xth)*log(1-s)+ApproxLogGenHarmSeries(1-s,N*xth,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,xth,s,a)
        y_aux = 1-((1-a)/(1-s))^(N*xth); % contributions to the cumulative p-distro starting below xc.
        y_exact = y_aux/(y_aux + (a-s)/(s*(1-a))*((1-a)/(1-s))^(N*xth)*((1/(1-s))^(N-N*xth)-1)   ); % contributions to the cumulative p-distro starting above xc.
        y_approx = 1/(1+ exp( N*xth*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,xth,s,a)
% Nc((Nc==0)|(Nc==N))=nan;
% Nth = N*xth;
% 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_approx = y_approx*(numel(Nc)-1)/sum(y_approx,'omitnan');
% 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)));
% y_exact = y_exact*(numel(Nc)-1)/sum(y_exact,'omitnan');
% end

% Function computing the mean extinction and fixation rates for 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_ext,y_fix] = ExtFixRates4FixXthSys(N,xth,s,a)
    [~,y1_exact] = FixProbAtThr4AsymDensDepSys(N,xth,s,a); % compute fixation probability (mean fixation rate divided by the sum of mean fixation and extinction rates).
    [~,y2_exact,~,~] = MeanAbsorptionTimeAtThr4AsymDensDepSys(N,xth,s,a); % compute mean coexistence time (inverse of the sum of mean fixation and extinction rates).
    
    y_ext = (1-y1_exact)./y2_exact; % overall mean extinction rate as the proportion of mean extinction rate with respect to the mean overall rate, multiplied by the mean overall rate.
    y_fix = y1_exact./y2_exact; % overall mean extinction rate as the proportion of mean fixation rate with respect to the mean overall rate, multiplied by the mean overall rate.
end

%% Function for the N-PDMP marginalisation
function [Y_approx,Y_exact,NPDMP] = NPDMPmarginalization(fhandle,xth,s,a,K,normalize)
    %fhandle=@funname;
    % x is a vector of the x-points evaluated that is used only when we marginalize de QSCD.
    % normalize is a logical value that either normalizes the result or doesn't when true or false, respectively.
    k1=K.values(1); k2=K.values(2); % Kminus and Kplus.
    Nuplus=K.nu.*(1-K.delta); % nu_down.
    Numinus=K.nu.*(1+K.delta); % nu_up.
    N_extrema_points_per_integer = 1e1; % the name explains it.
    N_points_per_integer = 1; % the name explains it.
    N = [((N_extrema_points_per_integer*k1):(N_extrema_points_per_integer*(k1+1)))/N_extrema_points_per_integer,((N_points_per_integer*(k1+2)):(N_points_per_integer*(k2-2)))/N_points_per_integer,((N_extrema_points_per_integer*(k2-1)):(N_extrema_points_per_integer*k2))/N_extrema_points_per_integer]; %((N_points_per_integer*k1):(N_points_per_integer*k2))/N_points_per_integer; % vector of N values.
    NPDMP = nan(2,numel(N),numel(Nuplus)); % initialize NPDMP distribution;
    % Check whether the chosen function input for NPDMPmarginalization uses the input 'x' and the output is a vector of size(x).
    if strcmp(func2str(fhandle),'FixAndCoexDiagram') % check whether NPDMPmarginalization uses the input 'x' and the output is a vector of size(x).
        Y_approx = nan(size(Nuplus)); Y_exact = Y_approx; % output dimensions;
        % Compute the function marginalized over N.
        f = nan(size(N)); g = f; f2=f; g2=f; % initialize the function evaluation.
        for N_idx = 1:numel(N) % for each (usually non-integer) total population in the [K-,K+] range.
            [g(N_idx),f(N_idx)] = FixProbAtThr4AsymDensDepSys(N(N_idx),xth,s,a); % compute the fixation probability
            [g2(N_idx),f2(N_idx)] = MeanAbsorptionTimeAtThr4AsymDensDepSys(N(N_idx),xth,s,a); % compute the mean coexistence time.
        end
        [Nmean,~]=NPDMPmarginalization(@(w,x,y,z)deal(w,0),[],s,a,K,true); % compute the mean total population for each environment.
        for env_idx=1:numel(Nuplus) % for each dynamic environment statistics nu and delta.
            % Compute the final function
            % extinction probability of resistant cells multiplied by non-coexistence probability.
            G=(1-f).*(1-exp(-2*Nmean(env_idx)./f2));%(f2<2*Nmean(env_idx));
            % fixation probability of resistant cells multiplied by non-coexistence probability.
            F=f.*(1-exp(-2*Nmean(env_idx)./f2));%(f2<2*Nmean(env_idx));
            % Compute the NPDMP for this environment.
            nuplus = Nuplus(env_idx); numinus = Numinus(env_idx); % nu rates.
            %Z = k2^numinus*k1^nuplus*gamma(nuplus+numinus)/( (k2-k1)^(nuplus+numinus-1)*gamma(nuplus)*gamma(numinus) );
            NPDMP_aux = (1-k1./N).^(numinus-1).*(k2./N-1).^(nuplus-1)./(N.^2); % use this N/k1 and N/k2 form so that it doesn't reach infinity.
            NPDMP(1,:,env_idx) = N;
            NPDMP(2,:,env_idx) = NPDMP_aux;
            % Compute the marginalization of the function F over the NPDMP avoiding numerical issues due to singularities at the bounds.
            y_k1 = zeros(2,1); y_k2 = y_k1; % auxiliar variables if there are singularities at the bounds.
            if numinus<1 % for the case that we have a singularity at kminus.
                k1plus1_idx = N==(k1+1); % select exactly the term that corresponds to N being kminus+1.
                idx_k1 = N<(k1+1); % select all terms with fractional population lower than kminus+1.
                N_k1 = N(idx_k1); % auxiliar vector for the fractional population values below kminus+1.
                NPDMP_aux_k1 = (1-k1./N_k1).^(numinus).*(k2./N_k1-1).^(nuplus-1)./(N_k1.^2); % NPDMP where the singularity has been integrated by parts.
                y_k1(1) = F(k1plus1_idx)*NPDMP_aux(k1plus1_idx)/numinus; % summand in the integration by parts.
                y_k1(1) = y_k1(1) - sum(((diff(F(idx_k1|k1plus1_idx))*N_extrema_points_per_integer - F(idx_k1).*((nuplus-1)./(k2-N_k1)+(nuplus+numinus)./N_k1))/numinus).*NPDMP_aux_k1); % total contribution of the integration by parts between kminus and kminus+1 (last point not included).
                y_k1(2) = G(k1plus1_idx)*NPDMP_aux(k1plus1_idx)/numinus; % summand in the integration by parts.
                y_k1(2) = y_k1(2) - sum(((diff(G(idx_k1|k1plus1_idx))*N_extrema_points_per_integer - G(idx_k1).*((nuplus-1)./(k2-N_k1)+(nuplus+numinus)./N_k1))/numinus).*NPDMP_aux_k1); % total contribution of the integration by parts between kminus and kminus+1 (last point not included).
            end
            if nuplus<1 % for the case that we have a singularity at kplus.
                k2minus1_idx = N==(k2-1); % select exactly the term that corresponds to N being kplus-1.
                idx_k2 = N>(k2-1); % select all terms with fractional population higher than kplus-1.
                N_k2 = N(idx_k2); % auxiliar vector for the fractional population values above kplus-1.
                NPDMP_aux_k2 = (1-k1./N_k2).^(numinus-1).*(k2./N_k2-1).^(nuplus)./(N_k2.^2); % NPDMP where the singularity has been integrated by parts.
                y_k2(1) = F(k2minus1_idx)*NPDMP_aux(k2minus1_idx)/nuplus; % summand in the integration by parts.
                y_k2(1) = y_k2(1) + sum(((flip(diff(flip(F(idx_k2|k2minus1_idx))))*N_extrema_points_per_integer + F(idx_k2).*((numinus-1)./(N_k2-k1)-(nuplus+numinus)./N_k2))/nuplus).*NPDMP_aux_k2); % total contribution of the integration by parts between kplus-1 and kplus (first point not included).
                y_k2(2) = G(k2minus1_idx)*NPDMP_aux(k2minus1_idx)/nuplus; % summand in the integration by parts.
                y_k2(2) = y_k2(2) + sum(((flip(diff(flip(G(idx_k2|k2minus1_idx))))*N_extrema_points_per_integer + G(idx_k2).*((numinus-1)./(N_k2-k1)-(nuplus+numinus)./N_k2))/nuplus).*NPDMP_aux_k2); % total contribution of the integration by parts between kplus-1 and kplus (first point not included).
            end
            idx_N = (N>=(k1+1))&(N<=(k2-1)); % for the inner domain, which has no singularities.
            Y_exact(env_idx) = sum(F(idx_N).*NPDMP_aux(idx_N))/N_points_per_integer + (y_k1(1) + y_k2(1))/N_extrema_points_per_integer;%sum(F(idx_N).*NPDMP_aux(idx_N)) + y_k1(1) + y_k2(1);
            Y_approx(env_idx) = sum(G(idx_N).*NPDMP_aux(idx_N))/N_points_per_integer + (y_k1(2) + y_k2(2))/N_extrema_points_per_integer;%sum(G(idx_N).*NPDMP_aux(idx_N)) + y_k1(2) + y_k2(2);
            if (numinus/nuplus)==0
                Y_exact(env_idx)=F(1); Y_approx(env_idx)=G(1);
            elseif (nuplus/numinus)==0
                Y_exact(env_idx)=F(end); Y_approx(env_idx)=G(end);
            end
        end
        if normalize
            [Normalization,~]=NPDMPmarginalization(@(w,x,y,z)deal(1,0),[],s,a,K,false);
            Y_exact = Y_exact./Normalization; Y_approx = Y_approx./Normalization;
            NPDMP(2,:,:) = squeeze(NPDMP(2,:,:))./Normalization;
        end
    elseif strcmp(func2str(fhandle),'AnnealedFixAndCoexDiagram')
        Y_approx = nan(size(Nuplus)); Y_exact = Y_approx; % output dimensions;
        % Compute the function marginalized over N.
        f = nan(size(N)); g = f; f2=f; g2=f; % initialize the function evaluation.
        for N_idx = 1:numel(N) % for each (usually non-integer) total population in the [K-,K+] range.
            [g(N_idx),f(N_idx)] = FixProbAtThr4AsymDensDepSys(N(N_idx),xth,s,a); % compute the fixation probability
            [g2(N_idx),f2(N_idx)] = MeanAbsorptionTimeAtThr4AsymDensDepSys(N(N_idx),xth,s,a); % compute the mean coexistence time.
        end
        [Nmean,~]=NPDMPmarginalization(@(w,x,y,z)deal(w,0),[],s,a,K,true); % compute the mean total population for each environment.
        for env_idx=1:numel(Nuplus) % for each dynamic environment statistics nu and delta.
            % Compute the final function
            % extinction probability of resistant cells multiplied by non-coexistence probability.
            G=(1-f)./f2.*(1-exp(-2*Nmean(env_idx)./f2));%(f2<2*Nmean(env_idx));
            % fixation probability of resistant cells multiplied by non-coexistence probability.
            F=f./f2.*(1-exp(-2*Nmean(env_idx)./f2));%(f2<2*Nmean(env_idx));
            % Compute the NPDMP for this environment.
            nuplus = Nuplus(env_idx); numinus = Numinus(env_idx); % nu rates.
            %Z = k2^numinus*k1^nuplus*gamma(nuplus+numinus)/( (k2-k1)^(nuplus+numinus-1)*gamma(nuplus)*gamma(numinus) );
            NPDMP_aux = (1-k1./N).^(numinus-1).*(k2./N-1).^(nuplus-1)./(N.^2); % use this N/k1 and N/k2 form so that it doesn't reach infinity.
            NPDMP(1,:,env_idx) = N;
            NPDMP(2,:,env_idx) = NPDMP_aux;
            % Compute the marginalization of the function F over the NPDMP avoiding numerical issues due to singularities at the bounds.
            y_k1 = zeros(2,1); y_k2 = y_k1; % auxiliar variables if there are singularities at the bounds.
            if numinus<1 % for the case that we have a singularity at kminus.
                k1plus1_idx = N==(k1+1); % select exactly the term that corresponds to N being kminus+1.
                idx_k1 = N<(k1+1); % select all terms with fractional population lower than kminus+1.
                N_k1 = N(idx_k1); % auxiliar vector for the fractional population values below kminus+1.
                NPDMP_aux_k1 = (1-k1./N_k1).^(numinus).*(k2./N_k1-1).^(nuplus-1)./(N_k1.^2); % NPDMP where the singularity has been integrated by parts.
                y_k1(1) = F(k1plus1_idx)*NPDMP_aux(k1plus1_idx)/numinus; % summand in the integration by parts.
                y_k1(1) = y_k1(1) - sum(((diff(F(idx_k1|k1plus1_idx))*N_extrema_points_per_integer - F(idx_k1).*((nuplus-1)./(k2-N_k1)+(nuplus+numinus)./N_k1))/numinus).*NPDMP_aux_k1); % total contribution of the integration by parts between kminus and kminus+1 (last point not included).
                y_k1(2) = G(k1plus1_idx)*NPDMP_aux(k1plus1_idx)/numinus; % summand in the integration by parts.
                y_k1(2) = y_k1(2) - sum(((diff(G(idx_k1|k1plus1_idx))*N_extrema_points_per_integer - G(idx_k1).*((nuplus-1)./(k2-N_k1)+(nuplus+numinus)./N_k1))/numinus).*NPDMP_aux_k1); % total contribution of the integration by parts between kminus and kminus+1 (last point not included).
            end
            if nuplus<1 % for the case that we have a singularity at kplus.
                k2minus1_idx = N==(k2-1); % select exactly the term that corresponds to N being kplus-1.
                idx_k2 = N>(k2-1); % select all terms with fractional population higher than kplus-1.
                N_k2 = N(idx_k2); % auxiliar vector for the fractional population values above kplus-1.
                NPDMP_aux_k2 = (1-k1./N_k2).^(numinus-1).*(k2./N_k2-1).^(nuplus)./(N_k2.^2); % NPDMP where the singularity has been integrated by parts.
                y_k2(1) = F(k2minus1_idx)*NPDMP_aux(k2minus1_idx)/nuplus; % summand in the integration by parts.
                y_k2(1) = y_k2(1) + sum(((flip(diff(flip(F(idx_k2|k2minus1_idx))))*N_extrema_points_per_integer + F(idx_k2).*((numinus-1)./(N_k2-k1)-(nuplus+numinus)./N_k2))/nuplus).*NPDMP_aux_k2); % total contribution of the integration by parts between kplus-1 and kplus (first point not included).
                y_k2(2) = G(k2minus1_idx)*NPDMP_aux(k2minus1_idx)/nuplus; % summand in the integration by parts.
                y_k2(2) = y_k2(2) + sum(((flip(diff(flip(G(idx_k2|k2minus1_idx))))*N_extrema_points_per_integer + G(idx_k2).*((numinus-1)./(N_k2-k1)-(nuplus+numinus)./N_k2))/nuplus).*NPDMP_aux_k2); % total contribution of the integration by parts between kplus-1 and kplus (first point not included).
            end
            idx_N = (N>=(k1+1))&(N<=(k2-1)); % for the inner domain, which has no singularities.
            Y_exact(env_idx) = sum(F(idx_N).*NPDMP_aux(idx_N))/N_points_per_integer + (y_k1(1) + y_k2(1))/N_extrema_points_per_integer;%sum(F(idx_N).*NPDMP_aux(idx_N)) + y_k1(1) + y_k2(1);
            Y_approx(env_idx) = sum(G(idx_N).*NPDMP_aux(idx_N))/N_points_per_integer + (y_k1(2) + y_k2(2))/N_extrema_points_per_integer;%sum(G(idx_N).*NPDMP_aux(idx_N)) + y_k1(2) + y_k2(2);
            if (numinus/nuplus)==0
                Y_exact(env_idx)=F(1); Y_approx(env_idx)=G(1);
            elseif (nuplus/numinus)==0
                Y_exact(env_idx)=F(end); Y_approx(env_idx)=G(end);
            end
        end
        if normalize
            [Normalization,~]=NPDMPmarginalization(@(w,x,y,z)deal(1,0),[],s,a,K,false);
            Y_exact = Y_exact./Normalization; Y_approx = Y_approx./Normalization;
            NPDMP(2,:,:) = squeeze(NPDMP(2,:,:))./Normalization;
        end
    else % NPDMPmarginalization ignores 'x' and gives a scalar output.
        Y_approx = nan(size(Nuplus)); Y_exact = Y_approx; % output dimensions;
        % Compute the function marginalized over N.
        F = nan(size(N)); G = F; % initialize the function evaluation.
        for N_idx = 1:numel(N) % for each (usually non-integer) total population in the [K-,K+] range.
            [G(N_idx),F(N_idx)] = fhandle(N(N_idx),xth,s,a); % compute the function.
        end
        for env_idx=1:numel(Nuplus) % for each dynamic environment statistics nu and delta.
            % Compute the NPDMP for this environment.
            nuplus = Nuplus(env_idx); numinus = Numinus(env_idx); % nu rates.
            %Z = k2^numinus*k1^nuplus*gamma(nuplus+numinus)/( (k2-k1)^(nuplus+numinus-1)*gamma(nuplus)*gamma(numinus) );
            NPDMP_aux = (1-k1./N).^(numinus-1).*(k2./N-1).^(nuplus-1)./(N.^2); % use this N/k1 and N/k2 form so that it doesn't reach infinity.
            NPDMP(1,:,env_idx) = N;
            NPDMP(2,:,env_idx) = NPDMP_aux;
            % Compute the marginalization of the function F over the NPDMP avoiding numerical issues due to singularities at the bounds.
            y_k1 = zeros(2,1); y_k2 = y_k1; % auxiliar variables if there are singularities at the bounds.
            if numinus<1 % for the case that we have a singularity at kminus.
                k1plus1_idx = N==(k1+1); % select exactly the term that corresponds to N being kminus+1.
                idx_k1 = N<(k1+1); % select all terms with fractional population lower than kminus+1.
                N_k1 = N(idx_k1); % auxiliar vector for the fractional population values below kminus+1.
                NPDMP_aux_k1 = (1-k1./N_k1).^(numinus).*(k2./N_k1-1).^(nuplus-1)./(N_k1.^2); % NPDMP where the singularity has been integrated by parts.
                y_k1(1) = F(k1plus1_idx)*NPDMP_aux(k1plus1_idx)/numinus; % summand in the integration by parts.
                y_k1(1) = y_k1(1) - sum(((diff(F(idx_k1|k1plus1_idx))*N_extrema_points_per_integer - F(idx_k1).*((nuplus-1)./(k2-N_k1)+(nuplus+numinus)./N_k1))/numinus).*NPDMP_aux_k1); % total contribution of the integration by parts between kminus and kminus+1 (last point not included).
                y_k1(2) = G(k1plus1_idx)*NPDMP_aux(k1plus1_idx)/numinus; % summand in the integration by parts.
                y_k1(2) = y_k1(2) - sum(((diff(G(idx_k1|k1plus1_idx))*N_extrema_points_per_integer - G(idx_k1).*((nuplus-1)./(k2-N_k1)+(nuplus+numinus)./N_k1))/numinus).*NPDMP_aux_k1); % total contribution of the integration by parts between kminus and kminus+1 (last point not included).
            end
            if nuplus<1 % for the case that we have a singularity at kplus.
                k2minus1_idx = N==(k2-1); % select exactly the term that corresponds to N being kplus-1.
                idx_k2 = N>(k2-1); % select all terms with fractional population higher than kplus-1.
                N_k2 = N(idx_k2); % auxiliar vector for the fractional population values above kplus-1.
                NPDMP_aux_k2 = (1-k1./N_k2).^(numinus-1).*(k2./N_k2-1).^(nuplus)./(N_k2.^2); % NPDMP where the singularity has been integrated by parts.
                y_k2(1) = F(k2minus1_idx)*NPDMP_aux(k2minus1_idx)/nuplus; % summand in the integration by parts.
                y_k2(1) = y_k2(1) + sum(((flip(diff(flip(F(idx_k2|k2minus1_idx))))*N_extrema_points_per_integer + F(idx_k2).*((numinus-1)./(N_k2-k1)-(nuplus+numinus)./N_k2))/nuplus).*NPDMP_aux_k2); % total contribution of the integration by parts between kplus-1 and kplus (first point not included).
                y_k2(2) = G(k2minus1_idx)*NPDMP_aux(k2minus1_idx)/nuplus; % summand in the integration by parts.
                y_k2(2) = y_k2(2) + sum(((flip(diff(flip(G(idx_k2|k2minus1_idx))))*N_extrema_points_per_integer + G(idx_k2).*((numinus-1)./(N_k2-k1)-(nuplus+numinus)./N_k2))/nuplus).*NPDMP_aux_k2); % total contribution of the integration by parts between kplus-1 and kplus (first point not included).
            end
            idx_N = (N>=(k1+1))&(N<=(k2-1)); % for the inner domain, which has no singularities.
            Y_exact(env_idx) = sum(F(idx_N).*NPDMP_aux(idx_N))/N_points_per_integer + (y_k1(1) + y_k2(1))/N_extrema_points_per_integer;%sum(F(idx_N).*NPDMP_aux(idx_N)) + y_k1(1) + y_k2(1);
            Y_approx(env_idx) = sum(G(idx_N).*NPDMP_aux(idx_N))/N_points_per_integer + (y_k1(2) + y_k2(2))/N_extrema_points_per_integer;%sum(G(idx_N).*NPDMP_aux(idx_N)) + y_k1(2) + y_k2(2);
            if (numinus/nuplus)==0
                Y_exact(env_idx)=F(1); Y_approx(env_idx)=G(1);
            elseif (nuplus/numinus)==0
                Y_exact(env_idx)=F(end); Y_approx(env_idx)=G(end);
            end
        end
        if normalize
            [Normalization,~]=NPDMPmarginalization(@(w,x,y,z)deal(1,0),[],s,a,K,false);
            Y_exact = Y_exact./Normalization; Y_approx = Y_approx./Normalization;
            NPDMP(2,:,:) = squeeze(NPDMP(2,:,:))./Normalization;
        end
    end
end