%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MATLAB Script to reproduce sea-level distributions from Coulon et al.
% (2023): calculation and plots of projections median and quantiles as well 
% as probability density functions at a given year for both prior and 
% posterior sea-level distributions.
%
% Author: Violaine Coulon (violaine.coulon@ulb.be)
% Date: October 2023
%
% Data Source:
% The data used in this script are available at 10.5281/zenodo.8398772
%
% MATLAB Version and Toolboxes:
% This script is designed to work with MATLAB R20XX or later and does not
% require any additional toolboxes beyond the standard MATLAB functions.
% Some functions used in the script might be available in earlier versions
% of MATLAB as well.
%
% Please make sure you have the necessary MATLAB version installed to run
% this script.
%
% Usage:
% 1. Ensure you have the necessary MATLAB version and toolboxes installed.
% 2. Download the paper and the required data from the provided sources.
% 3. Run the script 'Compute_Bayesian_Weight.m' to generate 'ENSEMBLE_SCORE.mat'
% 4. Run this script to generate the figure.
%
% Citation:
% If you use this script or the reproduced figure in your work, please
% consider citing the original paper as follows:
%
% Coulon, V., Klose, A. K., Kittel, C., Edwards, T., Turner, F., Winkelmann, R., and Pattyn, F.: 
% Disentangling the drivers of future Antarctic ice loss with a historically-calibrated ice-sheet model, 
% EGUsphere [preprint], https://doi.org/10.5194/egusphere-2023-1532, 2023.
%
% Disclaimer:
% This script is provided as-is without any warranties. The author is not
% responsible for any inaccuracies or errors in the reproduced figure.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;

fig=figure('Position', [100, 100, 1500, 800]);

subplot('Position', [0.05,0.55,0.26,0.4])
plot_SLC_mean_and_calibrated(2100)
ylim([-0.25 1.5])
legend({'SSP1-2.6','SSP5-8.5','SSP1-2.6 Bayesian Calibrated','SSP5-8.5 Bayesian Calibrated'},'Location','northwest','EdgeColor','none','FontSize',12)
subplot('Position', [0.38,0.55,0.26,0.4])
plot_SLC_mean_and_calibrated(2300)
legend off
ylim([-0.5 5])
subplot('Position', [0.71,0.55,0.26,0.4])
plot_SLC_mean_and_calibrated(3000)
legend off
ylim([-1.5 17])
yticks(0:2:18)

subplot('Position', [0.05,0.07,0.26,0.4])
plot_SLC_pdf_calibrated_vs_uncalibrated('SSP126',0,2100); hold on
plot_SLC_pdf_calibrated_vs_uncalibrated('SSP585',0,2100)
plot_SLC_pdf_calibrated_vs_uncalibrated('SSP126',1,2100); 
plot_SLC_pdf_calibrated_vs_uncalibrated('SSP585',1,2100)
legend off
subplot('Position', [0.38,0.07,0.26,0.4])
plot_SLC_pdf_calibrated_vs_uncalibrated('SSP126',0,2300); hold on
plot_SLC_pdf_calibrated_vs_uncalibrated('SSP585',0,2300)
plot_SLC_pdf_calibrated_vs_uncalibrated('SSP126',1,2300); 
plot_SLC_pdf_calibrated_vs_uncalibrated('SSP585',1,2300)
legend off
subplot('Position', [0.71,0.07,0.26,0.4])
plot_SLC_pdf_calibrated_vs_uncalibrated('SSP126',0,3000); hold on
plot_SLC_pdf_calibrated_vs_uncalibrated('SSP585',0,3000)
plot_SLC_pdf_calibrated_vs_uncalibrated('SSP126',1,3000); 
plot_SLC_pdf_calibrated_vs_uncalibrated('SSP585',1,3000)
xlim([-5 25])
legend off

set(gcf,'color','w');

function plot_SLC_pdf_calibrated_vs_uncalibrated(scenario,weighted,yr)

    load(['',scenario,'_ENSEMBLE_DATA.mat'],'SLC_ensemble','time')
    load ENSEMBLE_SCORE
    w1=w;
    
    if weighted==0
        w=zeros(1,size(SLC_ensemble,1))+1;
    end
    
    SLC_pdf_vs_time(SLC_ensemble,time,yr,scenario,w,weighted)

end

function SLC_pdf_vs_time(SLC_ensemble,time,time_snap,scenario,w,weighted)

    if isequal(scenario,'SSP585')
        cl=[178/255 34/255 34/255];
    else
        cl=[44/255 98/255 124/255];
    end
    if weighted==1
        typ='-';
        ll=2;
    else
        typ='--';
        ll=1;
    end
    for l=1:length(time_snap) % Go through target SLR level
        SLC=SLC_ensemble(:,find(time==time_snap(l)));
        [fi,xi] = ksdensity(SLC,'Weights',w);
        plot(xi,fi,'LineWidth',ll,'color',cl,'DisplayName',[num2str(time_snap(l)),' CE'],'LineStyle',char(typ)); hold on
    end
    ylabel('Probability Density')
    xlabel('Sea-Level Contribution (m)')

    set(gca,'FontSize',12)    
    legend('Location','NorthEast','FontSize',12,'EdgeColor','none')

end

function plot_SLC_mean_and_calibrated(end_time)

q=[0.05 0.5 0.95];

% BY SCENARIO

load ENSEMBLE_SCORE
w1=w;

%% ENSEMBLE MEAN

scenarios={'SSP126','SSP585'};

for weighted_mean=[0 1]
    for s=1:length(scenarios)

        scenario=char(scenarios(s));
        load(['Z:\Vio\PROTECT\RESULTS\ADDITIONAL_ENSEMBLES\',scenario,'_ENSEMBLE_DATA.mat'],'SLC_ensemble','time')
        if isequal(scenario,'SSP585')
            cl=[178/255 34/255 34/255];
        else
            cl=[44/255 98/255 124/255];
        end

        w=w1;
        if weighted_mean==0
            w=w.*0+1;
        end

        plot(time,time*0,'color',[0.5 0.5 0.5],'LineStyle',':','HandleVisibility','off'); hold on

        % Evaluate and plot quantile - 0.5-quantile = median value
        quart=zeros(length(q),size(SLC_ensemble,2));
        for ii=2:size(SLC_ensemble,2)
            quart(:,ii)=quantile_eval(SLC_ensemble(:,ii),w',q);
        end
        curve1=movmean(quart(1,:),5);
        curve2=movmean(quart(2,:),5);
        curve3=movmean(quart(3,:),5);

        time2 = [time, fliplr(time)];
        inBetween1 = [curve1, fliplr(curve3)];
        if weighted_mean==0
            h1=fill(time2, inBetween1,cl,'facealpha',.2,'EdgeColor',[0.8 0.8 0.8],'EdgeAlpha',.2,'HandleVisibility','off'); hold on
            h3=plot(time,curve2,'LineWidth',1,'Color',cl,'DisplayName',scenario,'LineStyle','--');
        else
            h1=fill(time2, inBetween1,cl,'facealpha',.4,'EdgeColor',[0.8 0.8 0.8],'EdgeAlpha',.4,'HandleVisibility','off'); hold on
            h3=plot(time,curve2,'LineWidth',2,'Color',cl,'DisplayName',[scenario,' Bayesian Calibrated']);
        end

    end
end

xlabel('Year (CE)')
ylabel('Sea-Level Contribution (m)')

legend('Location','NorthWest','FontSize',12,'EdgeColor','none','color','none')
set(gca,'FontSize',12)

xlim([2015 end_time])

if end_time==3000
    ylim([-2 16])
elseif end_time==2300
    ylim([-0.5 5])
end

end

% Evaluation of the alpha-quantile of the data set Y
% weight is a vector a weight assigned to the data
% alpha can be a vector

function quant = quantile_eval(Y,weight,alpha)

N = length(Y);

[Y_sorted,I] = sort(Y); % sort data in ascending order
weight = weight(I);
occ = countmember(Y_sorted,Y_sorted);
[Y_unique,ia] = unique(Y_sorted,'last');
pmasstilde = occ(ia)/N;
weighttilde = weight(ia);
Ftilde = cumsum(pmasstilde.*weighttilde);
Ftilde(end) = 1;
if(alpha>=Ftilde(end))
    quant = Y_unique(end);
    return
elseif(alpha<Ftilde(1))  
    quant = -inf;
    return
end 

quant = zeros(length(alpha),1);
% for ii=1:length(alpha)
%     ind = find(Ftilde>=alpha(ii),1,'first');
%     quant(ii) = Y_unique(ind);
% end
for ii=1:length(alpha)
  if(alpha(ii)>=Ftilde(end))
    quant(ii) = Y_unique(end);
  elseif(alpha(ii)<Ftilde(1))
    quant(ii) = -inf;
  else
    ind = find(Ftilde>=alpha(ii),1,'first');
    quant(ii) = Y_unique(ind);
  end
end

end

function COUNT = countmember(A,B)
% COUNTMEMBER - count members
%
%   COUNT = COUNTMEMBER(A,B) counts the number of times the elements of array A are
%   present in array B, so that C(k) equals the number of occurences of
%   A(k) in B. A may contain non-unique elements. C will have the same size as A.
%   A and B should be of the same type, and can be cell array of strings.
%

narginchk(2,2) ;
if ~isequal(class(A), class(B))
    error('Both inputs should be of the same class.') ;
end
if isempty(A) || isempty(B)
    % nothing to do
    COUNT = zeros(size(A)) ;
else    
    % which elements are unique in A,
    % also store the position to re-order later on
    [AUnique, ~, posJ] = unique(A(:)) ;
    % assign each element in B a number corresponding to the element of A
    [~, loc] = ismember(B, AUnique) ;
    % count these numbers
    N = histcounts(loc(:), 1:length(AUnique)+1) ;
    % re-order according to A, and reshape
    COUNT = reshape(N(posJ),size(A)) ;
end

end


