%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Figure 6a. TD/BU Driver Anomalies
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% NOTE: Must FIRST run Fig5a.m 
cd /Users/twrohr/Desktop/Matlab/Matlab/Grazing/CMIP6_Comp/NCC_Figs/Final_Figs
addpath(genpath('/Users/twrohr/Desktop/Matlab/Matlab/Utilities'))


% ===== Panels 2 & 3: TD−BU anomalies (metricIdx = 2), 
figure('Color','w','Units','normalized','Position',[0.04 0.10 1 1]);
W1 = 0.8; % width per panel 
H1 = .4; 
gh = .08; 
gw = .035; 
X0 = 0.07; % left margin 
Y0 = 0.1;

% ---- Move right panel further right (tweak) ----
dx  = 0.07;  % increase/decrease to taste
pos4 = [X0 + (W1-gw)./2+gw + dx,   Y0, (W1-gw)./2,  H1];

metricIdx = 2;
spread  = 0.42;                                  % half-width of model band around each driver
offsets = linspace(-0.5*spread, 0.5*spread, nmod);

pos1 = [X0,             Y0+H1+gh,         W1,  H1];                % Panel 1 (Mean -Absolute change)
pos3 = [X0 ,                  Y0+.05, (W1-gw)./2-.05,  H1+.4];                   % Panel 2 (TD–BU improving)
pos4 = [X0 + (W1-gw)./2+gw+.07,   Y0+.05, (W1-gw)./2-.05,  H1+.4];                   % Panel 3 (TD–BU worsening)


for p =  1:2
    % ----- Per-model means/std -----
    Tm = 100.*squeeze(Loss_Drivers(1:nmod, metricIdx, regTD(p), :, statMean));  % [nmod x 11]
    Bm = 100.*squeeze(Loss_Drivers(1:nmod, metricIdx, regBU(p), :, statMean));

    % Pair 11->6
    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) ];


    % Model anomaly = TDpair − BUpair
    Ymod_mean6 = TDpair_mean - BUpair_mean;

    % ----- Ensemble paired stacks (odd vs even components) -----
    Te = 100.*squeeze(Loss_Drivers(12, metricIdx, regTD(p), :, statMean));
    Be = 100.*squeeze(Loss_Drivers(12, metricIdx, regBU(p), :, statMean));

    lowSeg  = [ 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 = [ Te(2)-Be(2), Te(4)-Be(4), Te(6)-Be(6), Te(8)-Be(8), Te(10)-Be(10), 0 ];

    % ----- Axes for this panel -----
    posP = pos3; ttl = '\textbf{TD-BU (improving growth)}';
    if p==2, posP = pos4; ttl = '\textbf{TD-BU (worsening growth)}'; end

    % =========================
    % Axis A (LEFT): ensemble bars (normal sign)
    % =========================
    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];   % last is zero-height but keeps palette
    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,'XTickLabel',decomp_symbols3,'TickDir','out','FontSize',25);
    xtickangle(axL,0);

    yline(axL,0,'-','Color',[0 0 0],'LineWidth',0.9,'HandleVisibility','off');
    laxl = 100.*[-.25 .25];
    ylim(axL, laxl);

    %title(axL, ttl, 'Interpreter','latex');

    % Y-axis title #1 (left axis, this panel)
    ylabel(axL,{'\textbf{Top-down Anomaly}'},'Interpreter','latex');

    % =========================
    % Axis B (RIGHT): sign-flipped twin of Axis A (ticks/limits identical)
    % =========================
    axFlip = axes('Position', posP, 'Color','none', ...
                  'YAxisLocation','right', 'XAxisLocation','bottom', ...
                  'TickDir','out', 'Box','off');
    hold(axFlip,'on');

    xlim(axFlip, xlim(axL));
    ylim(axFlip, ylim(axL));

    % Hide x on this overlay axis so only axL shows x labels
    set(axFlip,'XTick',x6,'XTickLabel',[],'XColor','none','FontSize',25);

    yt = get(axL,'YTick');
    set(axFlip,'YTick',yt);
    set(axFlip,'YTickLabel', arrayfun(@(v) sprintf('%g', -v), yt, 'UniformOutput', false));

    % Y-axis title #2 (right sign-flipped axis, this panel)
    ylabel(axFlip,{'\textbf{Bottom-up Anomaly}'},'Interpreter','latex');

    % =========================
    % Axis C (MODELS): invisible overlay axis for stems + dots (no ticks/labels)
    % =========================
    axM = axes('Position', posP, 'Color','none', ...
               'YAxisLocation','right', 'XAxisLocation','bottom', ...
               'TickDir','out', 'Box','off');
    hold(axM,'on');

    % Match x; give models headroom like your original script
    xlim(axM, xlim(axL));
    ylim(axM, laxl.*2.5);

    % No ticks/labels (so we still have only the 4 “official” axes visible)
    set(axM,'XTick',[],'YTick',[],'XColor','none','YColor','none');

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

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

        % points
        scatter(axM, cx, Ymod_mean6(j,:), sz, c, 'filled', ...
            'MarkerEdgeColor','k', 'LineWidth',0.25);
    end

end


function y = pairsum(a,b)
    y = sum([a,b], 2, 'omitnan');        % row-wise sum ignoring NaNs
    y(isnan(a) & isnan(b)) = NaN;        % if both NaN -> NaN
end

