%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Supplementary Figure 2,12 - Model-specific biomass-weighted scatter plots
% One x-y variable pair; 11 models plus equal-model ensemble (4 x 3)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Currently Set up for Sup Fig 2; Adjust varaibles, names, bin number and axis limits to configure for Sup figure 12

clear; close all; clc;

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

% ======================================================================
% USER OPTIONS
% ======================================================================

% ---- Paths ----
dataDir = '/Users/twrohr/Desktop/data.nosync/SSP585/Abridged/regridded2';

% ---- Models ----
model_names = {'iHAMMOCC','CMOC','WOMBAT','OECOv2','MARBL','CanOE', ...
               'MEDUSA2.1','PISCESv2-IPSL','PISCESv2-CNRM', ...
               'BFM5.2','COBALTv2'};

mod_order  = [10,5,1,9,2,6,11,8,4,3,7];
modelsToUse = 1:11;

% ---- Single x-y variable pair ----
% x is plotted horizontally; y is plotted vertically.
xVar = struct('idx',6, ...
              'label','$\Delta_{\%} \mu  \; (\%)$');

yVar = struct('idx',1, ...
              'label','$\Delta_{\%} \Sigma P}  \; (\%)$'); % '$\Delta_{\%} \Sigma \textrm{Z}  \; (\%)$'

% ---- Optional TD-control / decoupled-biomass filter ----
% false = use all otherwise-valid biomass
% true  = retain cells where variables 1 and 6 have opposing signs
filterTDOnly = false;

% ---- Optional 1:1 line and six-region biomass fractions ----
% false = no 1:1 line; annotate the original four quadrant fractions
% true  = add y = x and split the upper-right and lower-left quadrants
%         at the 1:1 line, giving six biomass-weighted fractions
showOneToOne = false;

% ---- Optional biomass-weighted linear regression ----
% false = do not show regression
% true  = show the weighted least-squares line, approximate 95% confidence
%         interval, and weighted coefficient of determination in each panel
showRegression = true;

regLineColor    = [0 0 0];
regBandColor    = [0 0 0];
regLW           = 2.2;
regBandAlpha    = 0.006;
regTextFS       = 13;
regTextPosition = [0.025 0.52];  % normalized axes coordinates

% ---- Latitude filter ----
% Inclusive bounds in degrees north. Example Southern Ocean: [-80 -40].
minLat = -80;
maxLat =  80;

if minLat > maxLat
    error('minLat must be less than or equal to maxLat.');
end

% ---- Axes / binning ----
xL = [-100 100];
yL = [-100 100];
nbin = 150;

% Percentile trimming within each model; [0 100] means no trimming.
trimPct = [0 100];

% Remove a cell if either plotted value is below this threshold.
minPlotValue = -inf;

% ---- Density colours ----
% 'quadrant' = quadrant-specific RGB density map
% 'current'  = one conventional Crameri density map for every panel
colourMapMode = 'quadrant';

currentCmapName  = 'tofino';
currentCmapRange = [1 150];

clL = 1e-5;
clU = 0.015;

density_gamma = 1.2;
contrast_boost = 1.0;

showDensityColorbar = true;

% ---- Figure ----
figPosition = [0.03 0.03 0.58 0.92];
tileSpacing = 'compact';
figurePadding = 'loose';

% ======================================================================
% GLOBAL STYLE AND COLOURS
% ======================================================================

set(groot,'defaultAxesTickLabelInterpreter','latex', ...
          'defaultTextInterpreter','latex', ...
          'defaultLegendInterpreter','latex');

cmap1 = brewermap(13,'Oranges');
cmap2 = brewermap(15,'Greens');
cmap3 = brewermap(11,'Purples');

modelColors = [ ...
    cmap2([5,7,11,13,15],:); ...
    cmap1([5,7,11,13],:); ...
    cmap3([7,11],:); ...
    0 0 0 ...
];
modelColors = modelColors(1:11,:);

quadCols = makeQuadrantColours();
currentCmap = makeWhiteHeadCrameri( ...
    currentCmapName,currentCmapRange,50);

% ======================================================================
% LOAD AND CLEAN EACH MODEL ONCE
% ======================================================================

S = dir(fullfile(dataDir,'*.mat'));
S = S(~[S.isdir]);
[~,ord] = sort({S.name});
S = S(ord);

if numel(S) < max(mod_order(modelsToUse))
    error('Fewer model files were found than required by mod_order.');
end

panel = collectModelData( ...
    S,dataDir,mod_order,modelsToUse, ...
    xVar.idx,yVar.idx,trimPct,minPlotValue, ...
    filterTDOnly,minLat,maxLat);

% Equal-model ensemble: every valid model has total weight 1/nValid.
validModels = ~cellfun(@isempty,panel.xByModel);
nValidModels = sum(validModels);

if nValidModels == 0
    error('No models contain valid data after applying the filters.');
end

xEns = vertcat(panel.xByModel{validModels});
yEns = vertcat(panel.yByModel{validModels});
wEns = vertcat(panel.wByModel{validModels}) ./ nValidModels;

% ======================================================================
% 4 x 3 FIGURE: 11 MODELS + ENSEMBLE
% ======================================================================

f = figure('Units','normalized', ...
           'Position',figPosition, ...
           'Color','w');

tl = tiledlayout(f,4,3, ...
    'TileSpacing',tileSpacing, ...
    'Padding',figurePadding);

axH = gobjects(12,1);

for p = 1:12

    ax = nexttile(tl,p);
    axH(p) = ax;

    if p <= 11
        m = p;
        x = panel.xByModel{m};
        y = panel.yByModel{m};
        w = panel.wByModel{m};
        panelTitle = model_names{modelsToUse(m)};
        titleColor = modelColors(modelsToUse(m),:);
    else
        x = xEns;
        y = yEns;
        w = wEns;
        panelTitle = 'Ensemble';
        titleColor = [0 0 0];
    end

    plotWeightedPanel(ax,x,y,w,nbin,xL,yL, ...
        clL,clU,colourMapMode,currentCmap, ...
        density_gamma,contrast_boost,quadCols);

    hold(ax,'on');

    % Biomass-weighted linear regression with approximate 95% confidence
    % interval and weighted R^2 annotation.
    if showRegression
        addWeightedLinearRegression(ax,x,y,w,xL,yL, ...
            regLineColor,regBandColor,regLW,regBandAlpha, ...
            regTextFS,regTextPosition);
    end

    xline(ax,0,'k','LineWidth',1.5);
    yline(ax,0,'k','LineWidth',1.5);
    if showOneToOne
        addOneToOneLine(ax,xL,yL);
    end

    xlim(ax,xL);
    ylim(ax,yL);
    box(ax,'on');
    set(ax,'Layer','top','FontSize',12,'LineWidth',0.9);

    title(ax,panelTitle, ...
        'Color',titleColor, ...
        'Interpreter','latex', ...
        'FontWeight','bold', ...
        'FontSize',15);

    % Biomass-weighted fractions. With the 1:1 switch on, the same-sign
    % quadrants are each split at y = x, yielding six exhaustive regions.
    if showOneToOne
        addSixRegionPercW(ax,x,y,w,13);
    else
        addQuadPercW(ax,x,y,w,13);
    end
end

% Shared labels keep the 12 panels clean and aligned.
xlabel(tl,xVar.label,'Interpreter','latex','FontSize',19);
ylabel(tl,yVar.label,'Interpreter','latex','FontSize',19);

% Subplot letters a-l.
labs = arrayfun(@(k)sprintf('\\textbf{%s}',char('m'+k-1)), ...
                1:12,'UniformOutput',false);
placeSubplotLetters(axH,labs,[-0.03 1.015],17);

% One density key for the complete figure; no model contours or legend.
if showDensityColorbar
    drawnow;

    if strcmpi(colourMapMode,'quadrant')
        addQuadrantColorbarBesideAxes(f,axH(12),quadCols, ...
            clL,clU,density_gamma,contrast_boost,11);
    else
        cb = colorbar(axH(12));
        cb.FontSize = 11;
        cb.TickLabelInterpreter = 'latex';
        ylabel(cb,'Biomass-weighted frequency', ...
               'Interpreter','latex','FontSize',12);
    end
end

% ======================================================================
% HELPERS
% ======================================================================

function panel = collectModelData(S,dataDir,mod_order,modelsToUse, ...
                                  xIdx,yIdx,trimPct,minPlotValue, ...
                                  filterTDOnly,minLat,maxLat)

    nModels = numel(modelsToUse);
    xByModel = cell(nModels,1);
    yByModel = cell(nModels,1);
    wByModel = cell(nModels,1);

    for ii = 1:nModels

        m = modelsToUse(ii);
        fileNow = fullfile(dataDir,S(mod_order(m)).name);

        D = load(fileNow,'CoreVar_Inv_pc','Bio_Weights','Lat','Lon');

        Xmap = squeeze(D.CoreVar_Inv_pc(xIdx,:,:));
        Ymap = squeeze(D.CoreVar_Inv_pc(yIdx,:,:));
        Wmap = D.Bio_Weights;

        latMap = coordinateToMap(D.Lat,size(Wmap),'Lat');
        latitudeMask = isfinite(latMap) & ...
                       latMap >= minLat & latMap <= maxLat;

        Xmap(~latitudeMask) = NaN;
        Ymap(~latitudeMask) = NaN;
        Wmap(~latitudeMask) = NaN;

        if filterTDOnly
            biomassChange = squeeze(D.CoreVar_Inv_pc(1,:,:));
            growthChange  = squeeze(D.CoreVar_Inv_pc(6,:,:));

            tdMask = isfinite(biomassChange) & isfinite(growthChange) & ...
                     biomassChange .* growthChange < 0;

            Xmap(~tdMask) = NaN;
            Ymap(~tdMask) = NaN;
            Wmap(~tdMask) = NaN;
        end

        badCell = Xmap < minPlotValue | Ymap < minPlotValue;
        Xmap(badCell) = NaN;
        Ymap(badCell) = NaN;
        Wmap(badCell) = NaN;

        valid = isfinite(Xmap) & isfinite(Ymap) & ...
                isfinite(Wmap) & Wmap > 0;

        x = Xmap(valid);
        y = Ymap(valid);
        w = Wmap(valid);

        % Preserve relative biomass among cells, while setting every
        % model's total retained biomass weight to one.
        Wnorm = nansum(Wmap(:));
        if ~isfinite(Wnorm) || Wnorm <= 0
            continue
        end
        w = w ./ Wnorm;

        [x,y,w] = cleanAndTrimXYZ( ...
            x,y,w,trimPct,minPlotValue);

        xByModel{ii} = x;
        yByModel{ii} = y;
        wByModel{ii} = w;
    end

    panel.xByModel = xByModel;
    panel.yByModel = yByModel;
    panel.wByModel = wByModel;
end

function plotWeightedPanel(ax,x,y,w,nbin,xL,yL,clL,clU, ...
                           colourMapMode,currentCmap,density_gamma, ...
                           contrast_boost,quadCols)

    if isempty(x)
        cla(ax);
        text(ax,0.5,0.5,'No valid data', ...
            'Units','normalized','HorizontalAlignment','center', ...
            'Interpreter','latex','FontSize',12);
        xlim(ax,xL);
        ylim(ax,yL);
        return
    end

    switch lower(colourMapMode)
        case 'current'
            weightedBinscatter(ax,x,y,w,nbin,xL,yL);
            set(ax,'ColorScale','log');
            caxis(ax,[clL clU]);
            colormap(ax,currentCmap);

        case 'quadrant'
            plotQuadrantRGBDensity(ax,x,y,w,nbin,xL,yL, ...
                clL,clU,density_gamma,contrast_boost,quadCols);

        otherwise
            error('Unknown colourMapMode: %s.',colourMapMode);
    end
end

function coordMap = coordinateToMap(coord,targetSize,coordName)

    coord = squeeze(coord);

    if isequal(size(coord),targetSize)
        coordMap = coord;
    elseif isvector(coord) && numel(coord) == targetSize(1)
        coordMap = repmat(coord(:),1,targetSize(2));
    elseif isvector(coord) && numel(coord) == targetSize(2)
        coordMap = repmat(coord(:).',targetSize(1),1);
    elseif numel(coord) == prod(targetSize)
        coordMap = reshape(coord,targetSize);
    else
        error('%s dimensions [%s] do not match data dimensions [%s].', ...
              coordName,num2str(size(coord)),num2str(targetSize));
    end
end

function [x,y,w] = cleanAndTrimXYZ(x,y,w,trimPct,minPlotValue)

    x = x(:); y = y(:); w = w(:);
    good = isfinite(x) & isfinite(y) & isfinite(w) & w > 0 & ...
           x >= minPlotValue & y >= minPlotValue;
    x = x(good); y = y(good); w = w(good);

    if isempty(x)
        return
    end

    if trimPct(1) > 0 || trimPct(2) < 100
        px = prctile(x,trimPct);
        py = prctile(y,trimPct);
        keep = x >= px(1) & x <= px(2) & ...
               y >= py(1) & y <= py(2);
        x = x(keep); y = y(keep); w = w(keep);
    end
end

function cmap = makeWhiteHeadCrameri(cmapName,cmapRange,nWhite)

    baseN = cmapRange(2);
    crmap = crameri(cmapName,baseN);
    crmap = crmap(cmapRange(1):cmapRange(2),:);
    cmap = [ones(nWhite,3);crmap];

    for ii = 1:nWhite
        t = ii/nWhite;
        cmap(ii,:) = t.*cmap(nWhite+1,:) + (1-t).*[1 1 1];
    end
end

function weightedBinscatter(ax,x,y,w,nbin,xL,yL)

    [H,xedges,yedges] = weightedHist2D(x,y,w,nbin,xL,yL);
    H(H <= 0) = NaN;
    imagesc(ax,[xedges(1) xedges(end)], ...
               [yedges(1) yedges(end)],H.');
    set(ax,'YDir','normal');
    xlim(ax,xL); ylim(ax,yL);
end

function quadCols = makeQuadrantColours()

    white = [1.00 1.00 1.00];
    cool1 = [0.96 0.97 0.98];
    cool2 = [0.86 0.88 0.92];
    hot   = [0 0 0];

    blue = [0.70 0.83 0.98;0.43 0.65 0.90; ...
            0.22 0.45 0.76;0.08 0.24 0.55];
    green = [0.82 0.89 0.57;0.60 0.76 0.35; ...
             0.38 0.58 0.19;0.18 0.38 0.08];
    purple = [0.86 0.73 0.94;0.69 0.50 0.82; ...
              0.51 0.31 0.68;0.32 0.12 0.48];
    orange = [0.98 0.76 0.38;0.91 0.56 0.16; ...
              0.75 0.37 0.06;0.52 0.20 0.02];

    common = [white;cool1;cool2];
    quadCols.knots = [0 0.1 0.2 0.45 0.7 0.85 0.95 1];
    quadCols.BL_ramp = [common;blue;hot];
    quadCols.TR_ramp = [common;green;hot];
    quadCols.TL_ramp = [common;purple;hot];
    quadCols.BR_ramp = [common;orange;hot];
end

function plotQuadrantRGBDensity(ax,x,y,w,nbin,xL,yL, ...
                                clL,clU,density_gamma, ...
                                contrast_boost,quadCols)

    [H,xedges,yedges] = weightedHist2D(x,y,w,nbin,xL,yL);
    H(H <= 0) = NaN;
    Himg = H.';

    t = (log10(Himg)-log10(clL))./(log10(clU)-log10(clL));
    t = max(min(t,1),0).^density_gamma;
    t = min(1,contrast_boost.*t);

    xc = 0.5.*(xedges(1:end-1)+xedges(2:end));
    yc = 0.5.*(yedges(1:end-1)+yedges(2:end));
    [Xc,Yc] = meshgrid(xc,yc);

    rgb = ones(size(Himg,1),size(Himg,2),3);
    rgb = applyDensityRamp(rgb,t,Xc <= 0 & Yc <= 0, ...
                           quadCols.BL_ramp,quadCols.knots);
    rgb = applyDensityRamp(rgb,t,Xc > 0 & Yc > 0, ...
                           quadCols.TR_ramp,quadCols.knots);
    rgb = applyDensityRamp(rgb,t,Xc <= 0 & Yc > 0, ...
                           quadCols.TL_ramp,quadCols.knots);
    rgb = applyDensityRamp(rgb,t,Xc > 0 & Yc <= 0, ...
                           quadCols.BR_ramp,quadCols.knots);

    h = imagesc(ax,[xedges(1) xedges(end)], ...
                   [yedges(1) yedges(end)],rgb);
    set(h,'AlphaData',double(isfinite(Himg)));
    set(ax,'YDir','normal','Color',[1 1 1]);
    xlim(ax,xL); ylim(ax,yL);
end

function rgb = applyDensityRamp(rgb,t,qmask,ramp,knots)

    good = qmask & isfinite(t);
    if ~any(good(:)); return; end

    tt = t(good);
    for k = 1:3
        chan = rgb(:,:,k);
        chan(good) = interp1(knots,ramp(:,k),tt,'linear','extrap');
        rgb(:,:,k) = chan;
    end
    rgb = max(min(rgb,1),0);
end

function [H,xedges,yedges] = weightedHist2D(x,y,w,nbin,xL,yL)

    x = x(:); y = y(:); w = w(:);
    good = isfinite(x) & isfinite(y) & isfinite(w) & w > 0;
    x = x(good); y = y(good); w = w(good);

    xedges = linspace(xL(1),xL(2),nbin+1);
    yedges = linspace(yL(1),yL(2),nbin+1);

    ix = discretize(x,xedges);
    iy = discretize(y,yedges);
    keep = ~isnan(ix) & ~isnan(iy);

    if ~any(keep)
        H = zeros(nbin,nbin);
    else
        H = accumarray([ix(keep),iy(keep)],w(keep), ...
                       [nbin,nbin],@sum,0);
    end
end

function addQuadPercW(ax,x,y,w,fs)

    x = x(:); y = y(:); w = w(:);
    good = isfinite(x) & isfinite(y) & isfinite(w) & w > 0;
    x = x(good); y = y(good); w = w(good);

    W = max(nansum(w),1e-12);
    qTL = 100*nansum(w(x <= 0 & y > 0))/W;
    qTR = 100*nansum(w(x > 0 & y > 0))/W;
    qBR = 100*nansum(w(x > 0 & y <= 0))/W;
    qBL = 100*nansum(w(x <= 0 & y <= 0))/W;

    fmt = '\\textbf{%d\\%%}';
    m = 0.018;

    text(ax,m,1-m,sprintf(fmt,round(qTL)), ...
        'Units','normalized','HorizontalAlignment','left', ...
        'VerticalAlignment','top','FontSize',fs,'Interpreter','latex');
    text(ax,1-m,1-m,sprintf(fmt,round(qTR)), ...
        'Units','normalized','HorizontalAlignment','right', ...
        'VerticalAlignment','top','FontSize',fs,'Interpreter','latex');
    text(ax,1-m,m,sprintf(fmt,round(qBR)), ...
        'Units','normalized','HorizontalAlignment','right', ...
        'VerticalAlignment','bottom','FontSize',fs,'Interpreter','latex');
    text(ax,m,m,sprintf(fmt,round(qBL)), ...
        'Units','normalized','HorizontalAlignment','left', ...
        'VerticalAlignment','bottom','FontSize',fs,'Interpreter','latex');
end

function addOneToOneLine(ax,xL,yL)

    % Draw only the portion of y = x that lies inside both axis limits.
    lineL = [max(xL(1),yL(1)),min(xL(2),yL(2))];
    if lineL(1) <= lineL(2)
        plot(ax,lineL,lineL,'k--','LineWidth',1.25);
    end
end

function addSixRegionPercW(ax,x,y,w,fs)

    x = x(:); y = y(:); w = w(:);
    good = isfinite(x) & isfinite(y) & isfinite(w) & w > 0;
    x = x(good); y = y(good); w = w(good);

    W = max(nansum(w),1e-12);

    % Opposite-sign quadrants remain whole.
    qTL = 100*nansum(w(x <= 0 & y > 0))/W;
    qBR = 100*nansum(w(x > 0 & y <= 0))/W;

    % Split both same-sign quadrants at y = x. Equality is assigned to
    % the y <= x region so that all six fractions are exhaustive.
    qTR_above = 100*nansum(w(x > 0 & y > 0 & y >  x))/W;
    qTR_below = 100*nansum(w(x > 0 & y > 0 & y <= x))/W;
    qBL_above = 100*nansum(w(x <= 0 & y <= 0 & y >  x))/W;
    qBL_below = 100*nansum(w(x <= 0 & y <= 0 & y <= x))/W;

    fmt = '\\textbf{%d\\%%}';
    m = 0.018;

    % Whole opposite-sign quadrants.
    text(ax,m,1-m,sprintf(fmt,round(qTL)), ...
        'Units','normalized','HorizontalAlignment','left', ...
        'VerticalAlignment','top','FontSize',fs,'Interpreter','latex');
    text(ax,1-m,m,sprintf(fmt,round(qBR)), ...
        'Units','normalized','HorizontalAlignment','right', ...
        'VerticalAlignment','bottom','FontSize',fs,'Interpreter','latex');

    % Upper-right quadrant: above and below y = x.
    text(ax,0.62,1-m,sprintf(fmt,round(qTR_above)), ...
        'Units','normalized','HorizontalAlignment','center', ...
        'VerticalAlignment','top','FontSize',fs,'Interpreter','latex');
    text(ax,1-m,0.62,sprintf(fmt,round(qTR_below)), ...
        'Units','normalized','HorizontalAlignment','right', ...
        'VerticalAlignment','middle','FontSize',fs,'Interpreter','latex');

    % Lower-left quadrant: above and below y = x.
    text(ax,m,0.38,sprintf(fmt,round(qBL_above)), ...
        'Units','normalized','HorizontalAlignment','left', ...
        'VerticalAlignment','middle','FontSize',fs,'Interpreter','latex');
    text(ax,0.38,m,sprintf(fmt,round(qBL_below)), ...
        'Units','normalized','HorizontalAlignment','center', ...
        'VerticalAlignment','bottom','FontSize',fs,'Interpreter','latex');
end

function stats = addWeightedLinearRegression(ax,x,y,w,xL,yL, ...
                                             lineColor,bandColor, ...
                                             lineWidth,bandAlpha,fs, ...
                                             textPosition)
%ADDWEIGHTEDLINEARREGRESSION
% Add a biomass-weighted least-squares line, approximate 95% confidence
% interval for the mean fitted line, and weighted R^2 annotation.
%
% The fit is restricted to data inside the displayed axis limits.

    x = x(:);
    y = y(:);
    w = w(:);

    good = isfinite(x) & isfinite(y) & isfinite(w) & w > 0 & ...
           x >= xL(1) & x <= xL(2) & ...
           y >= yL(1) & y <= yL(2);

    x = x(good);
    y = y(good);
    w = w(good);

    stats = struct('slope',NaN,'intercept',NaN,'R2',NaN, ...
                   'slopeSE',NaN,'slopeCI95',NaN,'neff',NaN);

    if numel(x) < 5
        return
    end

    % Normalize biomass weights.
    w = w ./ nansum(w);

    % Effective sample size for the uncertainty estimate.
    neff = 1 ./ nansum(w.^2);
    p = 2;

    if ~isfinite(neff) || neff <= p + 1
        return
    end

    % Weighted least squares: y = beta(1) + beta(2)*x.
    X = [ones(size(x)),x];
    XtWX = X' * (X .* w);
    XtWy = X' * (y .* w);

    if rcond(XtWX) < eps
        return
    end

    beta = XtWX \ XtWy;

    yhat = X * beta;
    resid = y - yhat;

    % Weighted coefficient of determination.
    ybar = nansum(w .* y);
    SSE = nansum(w .* resid.^2);
    SST = nansum(w .* (y-ybar).^2);

    if ~isfinite(SST) || SST <= 0
        return
    end

    R2 = 1-SSE./SST;

    % Approximate residual variance, treating normalized biomass weights
    % as probability weights and using effective sample size.
    sigma2 = SSE .* neff ./ max(neff-p,1);
    covB = sigma2 .* inv(XtWX);

    slopeSE = sqrt(max(covB(2,2),0));

    if exist('tinv','file') == 2
        tcrit = tinv(0.975,max(neff-p,1));
    else
        tcrit = 1.96;
    end

    slopeCI95 = tcrit .* slopeSE;

    % Mean fitted response and its approximate 95% confidence interval.
    xg = linspace(xL(1),xL(2),300).';
    Xg = [ones(size(xg)),xg];
    yg = Xg * beta;

    varMean = sum((Xg*covB).*Xg,2);
    ci = tcrit .* sqrt(max(varMean,0));
    ylo = yg-ci;
    yhi = yg+ci;

    hold(ax,'on');

    if bandAlpha > 0
        fill(ax, ...
            [xg;flipud(xg)], ...
            [ylo;flipud(yhi)], ...
            bandColor, ...
            'EdgeColor','none', ...
            'FaceAlpha',bandAlpha, ...
            'HandleVisibility','off');
    end

    plot(ax,xg,ylo,'--', ...
        'Color',lineColor, ...
        'LineWidth',max(1,lineWidth*0.65), ...
        'HandleVisibility','off');

    plot(ax,xg,yhi,'--', ...
        'Color',lineColor, ...
        'LineWidth',max(1,lineWidth*0.65), ...
        'HandleVisibility','off');

    plot(ax,xg,yg,'-', ...
        'Color',lineColor, ...
        'LineWidth',lineWidth, ...
        'HandleVisibility','off');

    text(ax,textPosition(1),textPosition(2), ...
        sprintf('$R^2 = %.2f$',R2), ...
        'Units','normalized', ...
        'HorizontalAlignment','left', ...
        'VerticalAlignment','bottom', ...
        'FontSize',fs, ...
        'FontWeight','bold', ...
        'Interpreter','latex', ...
        'Color',lineColor);

    stats.slope = beta(2);
    stats.intercept = beta(1);
    stats.R2 = R2;
    stats.slopeSE = slopeSE;
    stats.slopeCI95 = slopeCI95;
    stats.neff = neff;
end

function cbAx = addQuadrantColorbarBesideAxes(parentFig,targetAx, ...
                                               quadCols,clL,clU, ...
                                               density_gamma, ...
                                               contrast_boost,fs)

    oldUnits = targetAx.Units;
    targetAx.Units = 'normalized';
    p = targetAx.Position;
    targetAx.Units = oldUnits;

    pos = [p(1)+p(3)+0.006,p(2),0.014,p(4)];
    n = 256;
    logVals = linspace(log10(clL),log10(clU),n).';
    t = (logVals-log10(clL))./(log10(clU)-log10(clL));
    t = min(1,contrast_boost.*(max(min(t,1),0).^density_gamma));

    ramps = {quadCols.BL_ramp,quadCols.TR_ramp, ...
             quadCols.TL_ramp,quadCols.BR_ramp};
    rgbBar = ones(n,4,3);

    for q = 1:4
        for k = 1:3
            rgbBar(:,q,k) = interp1(quadCols.knots,ramps{q}(:,k), ...
                                    t,'linear','extrap');
        end
    end

    cbAx = axes('Parent',parentFig,'Units','normalized','Position',pos);
    imagesc(cbAx,1:4,logVals,max(min(rgbBar,1),0));

    tickExp = ceil(log10(clL)):floor(log10(clU));
    tickVals = 10.^tickExp;
    set(cbAx,'YDir','normal','YAxisLocation','right', ...
        'XTick',[],'YTick',log10(tickVals), ...
        'YTickLabel',arrayfun(@(k)sprintf('$10^{%d}$',k), ...
                             tickExp,'UniformOutput',false), ...
        'TickLabelInterpreter','latex','FontSize',fs, ...
        'Box','on','LineWidth',0.8,'Layer','top');
    xlim(cbAx,[0.5 4.5]);
    ylim(cbAx,[log10(clL) log10(clU)]);
    ylabel(cbAx,'Biomass-weighted frequency', ...
           'Interpreter','latex','FontSize',fs+1);
end

function placeSubplotLetters(axList,labels,offset,fs)

    for k = 1:numel(axList)
        text(axList(k),offset(1),offset(2),labels{k}, ...
            'Units','normalized','Interpreter','latex', ...
            'HorizontalAlignment','left', ...
            'VerticalAlignment','bottom', ...
            'FontSize',fs,'Clipping','off');
    end
end