%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Figure 6. Magnitude Share Anomaly
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% IMPORTNAT! -- MUST FUNT Fig5a First!

%% %%%%%%%%%%%%
% Panel 2: TD - BU (Improving, worsening) u 
% ===== Panels 2 & 3: TD−BU anomalies (metricIdx = 2), overlay right axis for models =====
figure('Color','w','Units','normalized','Position',[0.04 0.10 1 1]);

% Three panels left-to-right
W1  = 0.8;          % width per panel
H1 = .4;
gh = .08;
gw = .035;
X0 = 0.06;          % left margin
Y0 = 0.1;

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


metricIdx = 2;
spread  = [.42];                                % half-width of model band around each driver
offsets = linspace(-.5*spread, .5*spread, nmod); % fixed left->right order


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));
    Ts = 100.*squeeze(Loss_Drivers(1:nmod, metricIdx, regTD(p), :, statStd));
    Bs = 100.*squeeze(Loss_Drivers(1:nmod, metricIdx, regBU(p), :, statStd));

    % Pair 11->6
    TDpair_mean = [ sum2(Tm(:,1),Tm(:,2)) , sum2(Tm(:,3),Tm(:,4)) , ...
                    sum2(Tm(:,5),Tm(:,6)) , sum2(Tm(:,7),Tm(:,8)) , ...
                    sum2(Tm(:,9),Tm(:,10)), Tm(:,11) ];
    BUpair_mean = [ sum2(Bm(:,1),Bm(:,2)) , sum2(Bm(:,3),Bm(:,4)) , ...
                    sum2(Bm(:,5),Bm(:,6)) , sum2(Bm(:,7),Bm(:,8)) , ...
                    sum2(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 ];
    Yens_stack_low  = lowSeg;
    Yens_stack_high = highSeg;

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

    % Left axis: ensemble stacks
    axL = axes('Position', posP); hold(axL,'on'); box(axL,'off');
    B = bar(axL, x6, [Yens_stack_low(:) Yens_stack_high(:)], '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,:);

    % Axis cosmetics (left / ensemble)
    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);
    ylabel(axL,'Ensemble Magnitude Share Anomaly ($\%$)');
    
    % Right overlay axis: model points (stems from 0 + dots)
    axR = axes('Position', posP, 'Color','none', 'YAxisLocation','right', ...
               'XAxisLocation','bottom', 'TickDir','out', 'Box','off'); 
    hold(axR,'on');

    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

    xlim(axR, xlim(axL));
    set(axR,'XTick',x6,'XTickLabel',[],'YTick',100.*[-.6:.2:.6],'FontSize',25);
    ylim(axR, laxl.*2.5);             % give models extra headroom if you like
    ylabel(axR,'Inter-model Magnitude Share Anomaly ($\%$)');

    % share x-axis; show ticks only on left axis unless p==2
    xlim(axR, xlim(axL));
    set(axR,'XTick',x6,'XTickLabel',[]);
    ylim(axR, laxl.*2.5);

    % Tick visibility rules:
    if p == 1
        % Hide RIGHT axis ticks/label on panel 2 (improving)
        set(axR, 'YTick', [], 'YColor', 'none');
        ylabel(axR, '');
    else
        % Hide LEFT axis ticks/label on panel 3 (worsening)
        set(axL, 'YTick', [], 'YColor', 'none');
        ylabel(axL, '');
    end

end




%% Support Functions
% -------------------- Helpers ------------------------
function [q50,q25,q75] = local_wquant3(x,w)
% Robust weighted 50/25/75 percentiles
    ok = isfinite(x) & isfinite(w) & (w > 0);
    x  = x(ok); w = w(ok);
    if isempty(x) || sum(w)<=0
        q50=NaN; q25=NaN; q75=NaN; return
    end
    [xs,ix] = sort(x(:));
    ws = w(:); ws = ws(ix);
    cw = cumsum(ws); cw = cw./cw(end);
    [cu,iu] = unique(cw,'stable'); xsu = xs(iu);
    if numel(cu)==1, q50=xsu; q25=xsu; q75=xsu; return, end
    q25 = interp1(cu, xsu, 0.25, 'linear','extrap');
    q50 = interp1(cu, xsu, 0.50, 'linear','extrap');
    q75 = interp1(cu, xsu, 0.75, 'linear','extrap');
end

% at the very end of the file (after your script code)
function y = sum2(a,b)
    y = sum([a,b], 2, 'omitnan');        % row-wise sum ignoring NaNs
    y(isnan(a) & isnan(b)) = NaN;        % if both NaN -> NaN
end