%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MATLAB Script to reproduce parameter probability distributions from Coulon 
% et al. (2023): calculation and plots of prior and posterior parameter 
% probability distributions functions (approximated using kernel density 
% estimation) and density histograms resulting from the calibration. 
%
% 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.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load ENSEMBLE_SCORE.mat
load LHSensemble.mat
nbin=20;
titles={'CLIM_{atm}','CLIM_{ocn}','\gamma_{atm} (°C km^{-1})','d_{ice} (m)','K_{ice} (mm PDD^{-1})','K_{snow} (mm PDD^{-1})','M_{param}','\Gamma_{eff}','GCM'};

fig=figure('Position', [200, 400, 1500, 600]);

for i=1:size(LHS,2)

    if i<3
        bin=2;
    elseif i==7
        bin=5;
    elseif i==9
        bin=4;
    else
        bin=nbin;
    end
    [bin_center,N]=weighted_histogram(LHS(:,i), w, bin)   ;

    subplot(2,5,i)
    if i==8
        bar(bin_center,N,'FaceColor',[32 104 139]/255,'EdgeColor',[32 104 139]/255,'LineWidth',2.5,'DisplayName','Calibrated histogram'); hold on
        legend('EdgeColor','none','FontSize',12,'Position',[0.81 0.25 0.06 0.06])
    else
        bar(bin_center,N,'HandleVisibility','off','FaceColor',[32 104 139]/255,'EdgeColor',[32 104 139]/255,'LineWidth',2.5); hold on
    end
    [fi,xi] = ksdensity(LHS(:,i),'Weights',w);
    pdi = makedist('Uniform');
    pdfi = pdf(pdi,xi);
    plot(xi,pdfi,'LineWidth',2,'DisplayName','Prior','LineStyle',':');
    if i==3 || i==4 || i==5 || i==6 || i==8
        plot(xi,fi,'LineWidth',2,'DisplayName','Posterior');
    end
    xlabel(char(titles(i)),'FontWeight','bold')
    ylim([0 2.3])
    set(gca,'ytick',[])
    xlim([xi(1) xi(end)])

    if i==1
        xticks([0.25 0.75]);
        xticklabels({'MAR','RACMO'})
        xlim([-0.05 1.05])
    elseif i==2
        xticks([0.25 0.75]);
        xticklabels({'ISMIP6','Schmidtko'})
        xlim([-0.05 1.05])
    elseif i==7
        xticks(bin_center);
        xticklabels({'M_{PICO}','M_{plume}','M_{quad}','M_{JD20}','M_{JD20s}'})
        xlim([-0.05 1.05])
    elseif i==9
        xticks(bin_center);
        xticklabels({'MRI-ESM','UKESM','CESM2','IPSL'})
        xlim([-0.05 1.05])
    elseif i<7
        x=linspace(0,1,5);
        xticks(x);
        mmax=max(LHval(:,i));
        mmin=min(LHval(:,i));
        x=x*(mmax-mmin)+mmin;
        if i==3
            x=round(x,4);
        else
            x=round(x,1);
        end
        xticklabels({num2str(x(1)),num2str(x(2)),num2str(x(3)),num2str(x(4)),num2str(x(5))})
    end

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

end

function [bin_center,N]=weighted_histogram(data,w,bin_number)
    
    [N,edges] = histcounts(data,bin_number); % histogram without weight for edges and normalization
    % define bin center
    bin_center=zeros(length(edges)-1,1);
    for b=1:length(edges)-1
    bin_center(b)=(edges(b)+edges(b+1))/2;
    end

    N_noweight=N; % save N without weight -- used for normalization
    N=N*0;

    % calculate N per bin
    for i=1:length(data)
        for b=1:length(edges)-1
            if data(i)>=edges(b) && data(i) <=edges(b+1)
                N(b)=N(b)+w(i);
            end
        end
    end

    N=N./N_noweight; % normalization for pdf

end