% This script produces figure 2
load('resilience_estimate_matrixP_high_res.mat')
load('resilience_estimate_matrixDD_high_res.mat')
%% 
figure
t= tiledlayout(2,1);
nexttile
h=imagesc(flip(resilience_estimate_matrix_p,1));
colormap parula(3);
cb = colorbar;
cb.Ticks = [0.3,1,1.7];
cb.TickLabels = {'Resilient', 'Criteria 2', 'Criteria 1'};
cb.Ruler.TickLabelRotation=270;
xticklabels = 0:10:100;
xticks = linspace(1, size(resilience_estimate_matrix_p, 2), numel(xticklabels));
set(gca, 'XTick', xticks, 'XTickLabel', xticklabels)
yticklabels = 0:10:100;
yticks = linspace(1, size(resilience_estimate_matrix_p, 1), numel(yticklabels));
set(gca, 'YTick', yticks, 'YTickLabel', flipud(yticklabels(:)))
xlabel('# Sites Disturbed (D)')
ylabel('# Sites Occupied (H)')
ax = gca;
ax.TitleHorizontalAlignment = 'left';
title('a) Proportional Kicks')

nexttile
h=imagesc(flip(resilience_estimate_matrix_dd,1));
colormap parula(3);
cb = colorbar;
cb.Ticks = [0.3,1,1.7];
cb.TickLabels = {'Resilient', 'Criteria 2', 'Criteria 1'};
cb.Ruler.TickLabelRotation=270;
xticklabels = 0:10:100;
xticks = linspace(1, size(resilience_estimate_matrix_dd, 2), numel(xticklabels));
set(gca, 'XTick', xticks, 'XTickLabel', xticklabels)
yticklabels = 0:10:100;
yticks = linspace(1, size(resilience_estimate_matrix_dd, 1), numel(yticklabels));
set(gca, 'YTick', yticks, 'YTickLabel', flipud(yticklabels(:)))
xlabel('# Sites Disturbed (D)')
ylabel('# Sites Occupied (H)')
ax = gca;
ax.TitleHorizontalAlignment = 'left';
title('b) Density-Dependent Kicks')