%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Figure 5a. TD/BU Driver Anomalies
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% NOTE: Must FIRST run Fig4a.m and not close/clc any filed

cd /Users/twrohr/Desktop/Matlab/Matlab/Grazing/CMIP6_Comp/NCC_Figs/Final_Figs/Revised
addpath(genpath('/Users/twrohr/Desktop/Matlab/Matlab/Utilities'))

% ============================================================
% Four-panel layout
%
% Standard quadrant numbering:
%
%       Q2 | Q1
%       -------
%       Q3 | Q4
%
% Quadrant reorganisation:
%   Original Q2 -> new Q1
%   Original Q1 -> new Q3
%   Original Q4 -> new Q2
%   Original Q3 -> new Q4
%
% Grey shading:
%   Left panels:  area below y = 0
%   Right panels: area above y = 0
%
% Left axis:  ensemble bars, ±25%
% Right axis: individual models, ±75%
% ============================================================

figure( ...
    'Color','w', ...
    'Units','normalized', ...
    'Position',[0.04 0.05 0.92 0.88]);

% ----- Analysis settings -----
metricIdx = 2;

spread = 0.42;
nmod   = 11;

offsets = linspace( ...
    -0.5*spread, ...
     0.5*spread, ...
     nmod);

% Sign applied to the original TD-BU anomaly:
% +1 = TD-BU
% -1 = BU-TD
rowSigns  = [-1 1];
rowLabels = {'TD-BU','BU-TD'};

% ----- Grey shading colour -----
shadeColor = [0.93 0.93 0.93];

% ----- Annotation settings -----
annotationFontSize = 16;
titleFontSize      = 25;
arrowMarkerSize    = 90;

% ============================================================
% Panel titles and directional annotations
% Indexed by displayed position: {row,column}
% ============================================================

panelTitles = { ...
    '\textbf{Decoupled Biomass Increase}', ...
    '\textbf{Coupled Biomass Increase}'; ...
    '\textbf{Coupled Biomass Decrease}', ...
    '\textbf{Decoupled Biomass Decrease}'};

topAnnotation = { ...
    {'Buffers','Decoupled','Increase'}, ...
    {'Buffers','Coupled','Increase'}; ...
    {'Amplifies','Coupled','Decrease'}, ...
    {'Favors','Decoupled','Decrease'}};

bottomAnnotation = { ...
    {'Favors','Decoupled','Increase'}, ...
    {'Amplifies','Coupled','Increase'}; ...
    {'Buffers','Coupled','Decrease'}, ...
    {'Buffers','Decoupled','Decrease'}};

% ----- Four-panel layout -----
leftMargin   = 0.035;
rightMargin  = 0.035;
bottomMargin = 0.05;
topMargin    = 0.175;

horizontalGap = 0.08;
verticalGap   = 0.05;

panelWidth = ...
    (1-leftMargin-rightMargin-horizontalGap)/2;

panelHeight = ...
    (1-bottomMargin-topMargin-verticalGap)/2;

xPos = [ ...
    leftMargin, ...
    leftMargin + panelWidth + horizontalGap];

yPos = [ ...
    bottomMargin + panelHeight + verticalGap, ... % top row
    bottomMargin];                                % bottom row

% ============================================================
% Destination-to-source quadrant mapping
%
% Destination position:
%
%   New Q2 (top-left)     <- original Q4
%   New Q1 (top-right)    <- original Q2
%   New Q3 (bottom-left)  <- original Q1
%   New Q4 (bottom-right) <- original Q3
%
% Each destination is indexed by:
%   row = 1 top,    row = 2 bottom
%   p   = 1 left,   p   = 2 right
% ============================================================

sourceRow = [ ...
    2 1; ...
    1 2];

sourceP = [ ...
    2 1; ...
    2 1];

for row = 1:2

    for p = 1:2

        % Original quadrant supplying data to this position
        dataRow = sourceRow(row,p);
        dataP   = sourceP(row,p);

        signFactor = rowSigns(dataRow);

        % =====================================================
        % Per-model means
        % =====================================================
        Tm = 100.*squeeze( ...
            Loss_Drivers( ...
            1:nmod, ...
            metricIdx, ...
            regTD(dataP), ...
            :, ...
            statMean));

        Bm = 100.*squeeze( ...
            Loss_Drivers( ...
            1:nmod, ...
            metricIdx, ...
            regBU(dataP), ...
            :, ...
            statMean));

        % Pair 11 components into 6 drivers
        TDpair_mean = [ ...
            pairsum(Tm(:,1),Tm(:,2)), ...
            pairsum(Tm(:,3),Tm(:,4)), ...
            pairsum(Tm(:,5),Tm(:,6)), ...
            pairsum(Tm(:,7),Tm(:,8)), ...
            pairsum(Tm(:,9),Tm(:,10)), ...
            Tm(:,11)];

        BUpair_mean = [ ...
            pairsum(Bm(:,1),Bm(:,2)), ...
            pairsum(Bm(:,3),Bm(:,4)), ...
            pairsum(Bm(:,5),Bm(:,6)), ...
            pairsum(Bm(:,7),Bm(:,8)), ...
            pairsum(Bm(:,9),Bm(:,10)), ...
            Bm(:,11)];

        % Original anomaly is TD-BU.
        % signFactor = -1 flips it to BU-TD.
        Ymod_mean6 = signFactor .* ...
            (TDpair_mean-BUpair_mean);

        % =====================================================
        % Ensemble paired stacks
        % =====================================================
        Te = 100.*squeeze( ...
            Loss_Drivers( ...
            12, ...
            metricIdx, ...
            regTD(dataP), ...
            :, ...
            statMean));

        Be = 100.*squeeze( ...
            Loss_Drivers( ...
            12, ...
            metricIdx, ...
            regBU(dataP), ...
            :, ...
            statMean));

        lowSeg = signFactor .* [ ...
            Te(1)-Be(1), ...
            Te(3)-Be(3), ...
            Te(5)-Be(5), ...
            Te(7)-Be(7), ...
            Te(9)-Be(9), ...
            Te(11)-Be(11)];

        highSeg = signFactor .* [ ...
            Te(2)-Be(2), ...
            Te(4)-Be(4), ...
            Te(6)-Be(6), ...
            Te(8)-Be(8), ...
            Te(10)-Be(10), ...
            0];

        % =====================================================
        % Panel position
        % =====================================================
        posP = [ ...
            xPos(p), ...
            yPos(row), ...
            panelWidth, ...
            panelHeight];

        % =====================================================
        % Axis A (LEFT): ensemble bars
        % =====================================================
        axL = axes('Position',posP);

        hold(axL,'on');
        box(axL,'off');

        B = bar( ...
            axL, ...
            x6, ...
            [lowSeg(:),highSeg(:)], ...
            'stacked', ...
            'BarWidth',0.72, ...
            'EdgeColor',[0 0 0], ...
            'LineWidth',1.0);

        oddMap  = [1 3 5 7 9 11];
        evenMap = [2 4 6 8 10 11];

        B(1).FaceColor = 'flat';
        B(2).FaceColor = 'flat';

        B(1).CData = com_col(oddMap,:);
        B(2).CData = com_col(evenMap,:);

        xlim(axL,[0.5-spread,6.5+spread]);

        set(axL, ...
            'XTick',x6, ...
            'TickDir','out', ...
            'FontSize',22, ...
            'YDir','normal');

        % Show x-axis labels only on bottom row
        if row == 2

            set(axL, ...
                'XTickLabel',decomp_symbols3);

            xtickangle(axL,0);

        else

            set(axL,'XTickLabel',[]);

        end

        % Ensemble y-axis limits: ±25%
        laxl = 100.*[-0.25 0.25];

        ylim(axL,laxl);

        % =====================================================
        % Light-grey background shading
        %
        % Both left panels:  shade below y = 0
        % Both right panels: shade above y = 0
        % =====================================================
        xl = xlim(axL);

        if p == 1

            % Left-column panels: lower half
            shadeY = [ ...
                laxl(1), ...
                laxl(1), ...
                0, ...
                0];

        else

            % Right-column panels: upper half
            shadeY = [ ...
                0, ...
                0, ...
                laxl(2), ...
                laxl(2)];

        end

        hShade = patch( ...
            axL, ...
            [xl(1),xl(2),xl(2),xl(1)], ...
            shadeY, ...
            shadeColor, ...
            'EdgeColor','none', ...
            'HandleVisibility','off', ...
            'HitTest','off');

        % Keep shading behind ensemble bars
        uistack(hShade,'bottom');

        yline(axL,0,'-', ...
            'Color',[0 0 0], ...
            'LineWidth',0.9, ...
            'HandleVisibility','off');

        % Optional left-axis label
        if p == 1

            % ylabel(axL, ...
            %     '\textbf{Ensemble anomaly (\%)}', ...
            %     'Interpreter','latex');

        end

        % =====================================================
        % Axis B (RIGHT): individual models
        % =====================================================
        axR = axes( ...
            'Position',posP, ...
            'Color','none', ...
            'YAxisLocation','right', ...
            'XAxisLocation','bottom', ...
            'TickDir','out', ...
            'Box','off', ...
            'YDir','normal');

        hold(axR,'on');

        % Match horizontal limits of Axis A
        xlim(axR,xlim(axL));

        % Models have three times the ensemble y-range
        ylim(axR,laxl.*3);

        % Axis A supplies visible x-axis information
        set(axR, ...
            'XTick',x6, ...
            'XTickLabel',[], ...
            'XColor','none', ...
            'YColor','k', ...
            'FontSize',22);

        yline(axR,0,'-', ...
            'Color',[0 0 0], ...
            'LineWidth',0.9, ...
            'HandleVisibility','off');

        % =====================================================
        % Plot model stems and points on Axis B
        % =====================================================
        for j = 1:nmod

            c  = cmap_mod(j,:);
            cx = x6 + offsets(j);

            % Stems
            line(axR, ...
                [cx;cx], ...
                [zeros(1,numel(cx));Ymod_mean6(j,:)], ...
                'Color',c, ...
                'LineWidth',lw, ...
                'HitTest','off');

            % Points
            scatter(axR, ...
                cx, ...
                Ymod_mean6(j,:), ...
                sz, ...
                c, ...
                'filled', ...
                'MarkerEdgeColor','k', ...
                'LineWidth',0.25);

        end

        % =====================================================
        % Panel title
        % =====================================================
        title(axR, ...
            panelTitles{row,p}, ...
            'Interpreter','latex', ...
            'FontSize',titleFontSize, ...
            'FontWeight','bold');

        % =====================================================
        % Directional interpretation annotations
        % =====================================================
        xlR = xlim(axR);
        ylR = ylim(axR);

        xRange = diff(xlR);
        yRange = diff(ylR);

        % Place text slightly inside left boundary
        xText = xlR(1) + 0.008*xRange;

        % Upper and lower text positions
        yTextTop = ...
            ylR(2) - 0.035*yRange;

        yTextBottom = ...
            ylR(1) + 0.035*yRange;

        text(axR, ...
            xText, ...
            yTextTop, ...
            topAnnotation{row,p}, ...
            'Interpreter','latex', ...
            'HorizontalAlignment','left', ...
            'VerticalAlignment','top', ...
            'FontSize',annotationFontSize, ...
            'Color','k', ...
            'Clipping','on');

        text(axR, ...
            xText, ...
            yTextBottom, ...
            bottomAnnotation{row,p}, ...
            'Interpreter','latex', ...
            'HorizontalAlignment','left', ...
            'VerticalAlignment','bottom', ...
            'FontSize',annotationFontSize, ...
            'Color','k', ...
            'Clipping','on');

        % Upward arrowhead
        scatter(axR, ...
            xlR(1), ...
            ylR(2), ...
            arrowMarkerSize, ...
            'k', ...
            '^', ...
            'filled', ...
            'Clipping','off', ...
            'HandleVisibility','off');

        % Downward arrowhead
        scatter(axR, ...
            xlR(1), ...
            ylR(1), ...
            arrowMarkerSize, ...
            'k', ...
            'v', ...
            'filled', ...
            'Clipping','off', ...
            'HandleVisibility','off');

        % Optional right-axis label
        if p == 2

            % ylabel(axR, ...
            %     '\textbf{Individual-model anomaly (\%)}', ...
            %     'Interpreter','latex');

        end

        % Keep axis lines and ticks above plotted objects
        axL.Layer = 'top';
        axR.Layer = 'top';

    end

end

function y = pairsum(a,b)

    % Row-wise sum, ignoring individual NaNs
    y = sum([a,b],2,'omitnan');

    % Preserve NaN if both components are NaN
    y(isnan(a) & isnan(b)) = NaN;

end