%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Sup Figure 14 - Bi Variate Scatter plots for Six independently configurable x-y variable pairs; optional TD mask
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear; close all; clc;

cd /Users/twrohr/Desktop/Matlab/Matlab/Grazing/CMIP6_Comp/NCC_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'};

% File order in dataDir, mapped onto model_names order above
mod_order = [10,5,1,9,2,6,11,8,4,3,7];
modelsToUse = 1:11;   % all
modelsToUse = [8,11]; % Just fishmip

% ---- Optional TD-control / decoupled-biomass filter ----
% false = use all otherwise-valid biomass (current behaviour)
% true  = retain only grid cells where variables 1 and 6 have opposing
%         signs: Delta SigmaP .* Delta mu < 0
filterTDOnly = false;

% ---- Latitude filter ----
% Retain grid cells within these inclusive latitude bounds (degrees north).
% Use [-90 90] for the full globe; e.g., [-90 -30] for south of 30 S.
minLat = -80;
maxLat =  -40;

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

% ---- Six independent x-y variable pairs (shown row-wise in a 3x2 grid) ----
% Edit xIdx and yIdx independently for every panel. Each panel also has its
% own title, axis labels, and (when colourMapMode='current') density map.
panelVars = [ ...
    struct('xIdx',1, 'yIdx',50, ...
      'title','Panel 1', ...
      'xLabel','$\Delta_{\%} \textrm{Phytoplankton Biomass}\; (\%)$', ...
      'yLabel','$\Delta_{\%} \textrm{Trophic Transfer} \; (\%)$', ...
      'cmap','tofino','cmapRange',[1 150]); ...
   struct('xIdx',3, 'yIdx',5, ...
      'title','Panel 5', ...
      'xLabel','$\Delta_{\%} \textrm{Net Primary Production} \; (\%)$', ...
      'yLabel','$\Delta_{\%} \textrm{Trophic Transfer} \; (\%)$', ...
      'cmap','lisbon','cmapRange',[21 170]); ... 
    struct('xIdx',2, 'yIdx',50, ...
      'title','Panel 2', ...
      'xLabel','$\Delta_{\%} \textrm{Zooplankton Biomass}\; (\%)$', ...
      'yLabel','$\Delta_{\%} \textrm{Trophic Transfer} \; (\%)$', ...
      'cmap','tofino','cmapRange',[1 150]); ...
    struct('xIdx',8, 'yIdx',50, ...
      'title','Panel 4', ...
      'xLabel','$\Delta_{\%} \textrm{Gross Secondary Production} \; (\%)$', ...
      'yLabel','$\Delta_{\%} \textrm{Trophic Transfer} \; (\%)$', ...
      'cmap','berlin','cmapRange',[1 150]); ... 
    struct('xIdx',52, 'yIdx',50, ...
      'title','Panel 3', ...
      'xLabel','$\Delta_{\%} \textrm{Phyto + Zooplankton Biomass}\; (\%)$', ...
      'yLabel','$\Delta_{\%} \textrm{Trophic Transfer} \; (\%)$', ...
      'cmap','berlin','cmapRange',[1 150]); ...
    struct('xIdx',49, 'yIdx',50, ...
      'title','Panel 6', ...
      'xLabel','$\Delta_{\%} \textrm{Mesozooplankton Production} \; (\%)$', ...
      'yLabel','$\Delta_{\%} \textrm{Trophic Transfer} \; (\%)$', ...
      'cmap','lisbon','cmapRange',[21 170]) ...
    % struct('xIdx',51, 'yIdx',50, ...
    %   'title','Panel 6', ...
    %   'xLabel','$\Delta_{\%} \textrm{MQuadratic Mortality} \; (\%)$', ...
    %   'yLabel','$\Delta_{\%} \textrm{Trophic Transfer} \; (\%)$', ...
    %   'cmap','lisbon','cmapRange',[21 170]) ...
];

if numel(panelVars) ~= 6
    error('panelVars must contain exactly six x-y variable definitions.');
end

% ---- Axes / binning ----
xL = [-90 90];
yL = [-90 90];

nbin = 150;

% Percentile trimming applied within each model before combining.
% Use [0 100] for effectively no trimming.
trimPct = [0 100];

% Remove a grid cell from a subplot if either plotted variable is below
% this value. Applied separately to each model and x/y subplot pairing.
minPlotValue = -inf;%100;

% Colour limits for weighted bin mass
clL = 1e-5;
clU = 0.015;

% ---- Density colour-map mode ----
% 'current'  = retain the existing row-specific Crameri maps exactly
% 'quadrant' = use the quadrant-specific RGB density map from Figure 1e
colourMapMode = 'quadrant';

% These settings are used only when colourMapMode = 'quadrant'.
density_gamma = 1.2;
contrast_boost = 1.0;

% Quadrant RGB ramps: BL blue, TR green, TL purple, BR orange.
quadCols = makeQuadrantColours();

% ---- Per-model contour overlay ----
showModelContours = true;     % switch false to recover original no-contour look
contourMassLevel  = 0.9;     % contour enclosing densest 70% of model-weighted mass
contourNBins      = 150;
contourSmoothSig  = 4;      % smoothing in bin units; 0 = no smoothing
contourLW         = .7;
contourMinPoints  = 150;

% Optional legend for contours. Usually cleaner left off.
showContourLegend = true;

% ---- Biomass-weighted linear regression ----
showRegression = true;
regressionConfidence = 0.95;  % confidence interval for the mean fitted line
regressionLineColor = [0.10 0.10 0.10];
regressionBandColor = [0.25 0.25 0.25];
regressionBandAlpha = 0.18;
regressionLineWidth = 1.8;
regressionLabelXFraction = 0.72; % position along fitted x range, from 0 to 1
regressionFontSize = 13;

% ---- Figure ----
figPosition = [0.05 0.05 0.5 1];

% ======================================================================
% GLOBAL STYLE
% ======================================================================

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

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

tl = tiledlayout(f, 3, 2, ...
                 'TileSpacing','compact', ...
                 'Padding','loose');

% ======================================================================
% MODEL COLOURS FROM OTHER SCRIPTS
% ======================================================================

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,:);

% ======================================================================
% FILE LIST
% ======================================================================

S = dir(dataDir);
S = S(~[S.isdir]);

% MATLAB usually returns dir alphabetically, but make this explicit.
[~,ord] = sort({S.name});
S = S(ord);

% ======================================================================
% MAIN LOOP
% ======================================================================

nRows = 3;
nCols = 2;

axH = gobjects(nRows,nCols);

legendHandles = gobjects(numel(modelsToUse),1);

for r = 1:nRows
    for c = 1:nCols

        panelNumber = (r-1)*nCols + c;
        thisPanel = panelVars(panelNumber);

        % ---- Panel-specific colormap with white head ----
        cmap = makeWhiteHeadCrameri( ...
            thisPanel.cmap, thisPanel.cmapRange, 50);

        ax = nexttile(tl, (r-1)*nCols + c);
        axH(r,c) = ax;

        % ---- Gather model data for this x/y pair ----
        panel = collectPanelData( ...
            S, dataDir, mod_order, modelsToUse, ...
            thisPanel.xIdx, thisPanel.yIdx, trimPct, minPlotValue, ...
            filterTDOnly, minLat, maxLat);

        % ---- Weighted ensemble binscatter ----
        switch lower(colourMapMode)
            case 'current'
                % Original plotting configuration, unchanged.
                weightedBinscatter(ax, panel.xAll, panel.yAll, panel.wAll, ...
                                   nbin, xL, yL);

                set(ax,'ColorScale','log');
                caxis(ax, [clL clU]);
                colormap(ax, cmap);

            case 'quadrant'
                plotQuadrantRGBDensity(ax, panel.xAll, panel.yAll, panel.wAll, ...
                    nbin, xL, yL, clL, clU, density_gamma, ...
                    contrast_boost, quadCols);

            otherwise
                error('Unknown colourMapMode: %s. Use ''current'' or ''quadrant''.', ...
                      colourMapMode);
        end

        hold(ax,'on');

        % ---- Biomass-weighted linear regression and confidence band ----
        if showRegression
            addWeightedLinearRegression(ax, ...
                panel.xAll, panel.yAll, panel.wAll, xL, yL, ...
                regressionConfidence, regressionLineColor, ...
                regressionBandColor, regressionBandAlpha, ...
                regressionLineWidth, regressionLabelXFraction, ...
                regressionFontSize);
        end

        % ---- Optional model contours ----
        if showModelContours
            hTmp = plotModelContours(ax, panel, modelColors, ...
                                     xL, yL, contourNBins, ...
                                     contourMassLevel, contourSmoothSig, ...
                                     contourLW, contourMinPoints);

            if r == 1 && c == 1
                legendHandles = hTmp;
            end
        end

        % ---- Zero lines ----
        xline(ax,0,'k','LineWidth',2);
        yline(ax,0,'k','LineWidth',2);

        xlim(ax,xL);
        ylim(ax,yL);

        box(ax,'on');
        set(ax,'Layer','top');

        % ---- Titles / labels ----
        % title(ax, thisPanel.title, 'Interpreter','latex','FontSize',20);
        xlabel(ax, thisPanel.xLabel, 'FontSize',16);
        ylabel(ax, thisPanel.yLabel, 'FontSize',16);

        % ---- Quadrant percentages ----
        addQuadPercW(ax, panel.xAll, panel.yAll, panel.wAll);

        % ---- Colourbar only on right column ----
        if c == nCols && strcmpi(colourMapMode,'current')
            cb = colorbar(ax);
            cb.FontSize = 14;
            cb.TickLabelInterpreter = 'latex';
            ylabel(cb,'Biomass-weighted frequency', ...
                   'Interpreter','latex', ...
                   'FontSize',16);

            tmin = floor(log10(clL));
            tmax = ceil(log10(clU));
            cb.Ticks = 10.^(tmin:tmax);
            cb.TickLabels = arrayfun(@(k) sprintf('$10^{%d}$',k), ...
                                      tmin:tmax, ...
                                      'UniformOutput',false);
        elseif strcmpi(colourMapMode,'current')
            colorbar(ax,'off');
        end
    end
end

% The quadrant maps are RGB images, so MATLAB's ordinary colorbar cannot
% represent them. Add the same four-strip density key used in Figure 1e.
if strcmpi(colourMapMode,'quadrant')
    drawnow;
    for r = 1:nRows
        addQuadrantColorbarBesideAxes(f, axH(r,nCols), quadCols, ...
            clL, clU, density_gamma, contrast_boost, 14);
    end
end

% ---- Optional model contour legend: coloured text, two bottom rows ----
if showModelContours && showContourLegend

    goodModels = modelsToUse(isgraphics(legendHandles));
    nModels    = numel(goodModels);

    nRowsLeg = 2;
    nColsLeg = ceil(nModels/nRowsLeg);

    % Legend positioning in normalized figure coordinates
    xLeft  = 0.1;
    xRight = 0.94;
    yPos   = [0.030, 0.008];   % top and bottom legend rows
    boxH   = 0.022;

    colW = (xRight-xLeft)/nColsLeg;

    for ii = 1:nModels

        % Fill the first row, then the second row
        legRow = ceil(ii/nColsLeg);
        legCol = mod(ii-1,nColsLeg) + 1;

        m = goodModels(ii);

        xPos = xLeft + (legCol-1)*colW;

        annotation(f,'textbox', ...
            [xPos, yPos(legRow), colW, boxH], ...
            'String',model_names{m}, ...
            'Color',modelColors(m,:), ...
            'Interpreter','latex', ...
            'FontSize',15, ...
            'FontWeight','bold', ...
            'HorizontalAlignment','center', ...
            'VerticalAlignment','middle', ...
            'LineStyle','none', ...
            'FitBoxToText','off');
    end
end

% ---- Subplot letters a-f ----
axList = axH.';
axList = axList(:);

labs = arrayfun(@(k) sprintf('\\textbf{%s}', char('a'+k-1)), ...
                1:numel(axList), ...
                'UniformOutput',false);

%placeSubplotLetters(axList, labs, [-0.115 .99], 24);

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

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

    xAll = [];
    yAll = [];
    wAll = [];

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

    for ii = 1:numel(modelsToUse)

        m = modelsToUse(ii);
        modelIndex(ii) = m;

        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;

        % Convert Lat to a 2-D map matching the plotted data. This supports
        % either a full latitude grid or a 1-D latitude coordinate vector.
        latMap = coordinateToMap(D.Lat,size(Wmap),'Lat');

        % Apply the latitude filter before model normalization or any panel
        % calculations, so excluded cells cannot affect density, quadrant
        % percentages, or model-specific contours.
        latitudeMask = isfinite(latMap) & ...
                       latMap >= minLat & latMap <= maxLat;

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

        % Optional common TD-control mask. This is always defined from
        % variables 1 (Delta SigmaP) and 6 (Delta mu), regardless of which
        % variables are plotted in the current subplot. Opposing signs are
        % treated as decoupled / top-down-controlled biomass.
        if filterTDOnly
            % decoupling mask 
            biomassChange = squeeze(D.CoreVar_Inv_pc(1,:,:));
            growthChange  = squeeze(D.CoreVar_Inv_pc(6,:,:));

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

            % % Community Comp mask
            % SP_Change  = squeeze(D.CoreVar_Inv_pc(12,:,:));
            % tdMask = isfinite(SP_Change)  & ...
            %          SP_Change < 0;


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

        % Clean the loaded model arrays before any downstream calculation.
        % If either plotted variable is below the threshold, invalidate the
        % complete grid cell, including its biomass weight. Consequently it
        % cannot enter normalization, density, quadrant percentages, or the
        % model-specific contours.
        badCell = Xmap < minPlotValue | Ymap < minPlotValue;
        Xmap(badCell) = NaN;
        Ymap(badCell) = NaN;
        Wmap(badCell) = NaN;

        % Use every remaining grid cell with valid x/y/weight.
        M = isfinite(Xmap) & isfinite(Ymap) & ...
            isfinite(Wmap) & Wmap > 0;

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

        % Normalize after cleaning, so each model's remaining valid biomass
        % weight sums to one.
        Wnorm = nansum(Wmap(:));
        if ~isfinite(Wnorm) || Wnorm <= 0
            xByModel{ii} = [];
            yByModel{ii} = [];
            wByModel{ii} = [];
            continue
        end

        w = w ./ Wnorm;

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

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

        xAll = [xAll; x];
        yAll = [yAll; y];
        wAll = [wAll; w];
    end

    panel.xAll = xAll;
    panel.yAll = yAll;
    panel.wAll = wAll;

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

function coordMap = coordinateToMap(coord,targetSize,coordName)
%COORDINATETOMAP Convert a coordinate grid/vector to the data-map shape.

    coord = squeeze(coord);

    if isequal(size(coord),targetSize)
        coordMap = coord;

    elseif isvector(coord) && numel(coord) == targetSize(1)
        % Latitude-style vector varying down the rows.
        coordMap = repmat(coord(:),1,targetSize(2));

    elseif isvector(coord) && numel(coord) == targetSize(2)
        % Coordinate vector varying across the columns.
        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(:);

    % The threshold cleaner has already been applied to the full map arrays.
    % Retain this threshold check as a defensive safeguard.
    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, yedges, H.');
    set(ax,'YDir','normal');

    xlim(ax,xL);
    ylim(ax,yL);
end

function quadCols = makeQuadrantColours()

    white_rgb = [1.00, 1.00, 1.00];
    cool1_rgb = [0.96, 0.97, 0.98];
    cool2_rgb = [0.86, 0.88, 0.92];
    hot_rgb   = [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];

    commonRamp = [white_rgb; cool1_rgb; cool2_rgb];

    quadCols.knots  = [0 0.1 0.2 0.45 0.7 0.85 0.95 1];
    quadCols.BL_ramp = [commonRamp; blue;   hot_rgb];
    quadCols.TR_ramp = [commonRamp; green;  hot_rgb];
    quadCols.TL_ramp = [commonRamp; purple; hot_rgb];
    quadCols.BR_ramp = [commonRamp; orange; hot_rgb];
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);
    t = t .^ 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 cbAx = addQuadrantColorbarBesideAxes(parentFig, targetAx, quadCols, ...
                                               clL, clU, density_gamma, ...
                                               contrast_boost, fs)

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

    gap = 0.008;
    barWidth = 0.018;
    pos = [p(1)+p(3)+gap, p(2), barWidth, 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',[],'XTickLabel',[], ...
        'YTick',log10(tickVals), ...
        'YTickLabel',arrayfun(@(k)sprintf('$10^{%d}$',k),tickExp, ...
                             'UniformOutput',false), ...
        'TickLabelInterpreter','latex','FontSize',fs, ...
        'Box','on','LineWidth',1,'Layer','top');
    xlim(cbAx,[0.5 4.5]);
    ylim(cbAx,[log10(clL) log10(clU)]);
    ylabel(cbAx,'Biomass-weighted frequency', ...
           'Interpreter','latex','FontSize',16);
end

function hOut = plotModelContours(ax, panel, modelColors, xL, yL, nbin, ...
                                  massLevel, smoothSig, lineWidth, minPoints)

    hOut = gobjects(numel(panel.xByModel),1);

    hold(ax,'on');

    for ii = 1:numel(panel.xByModel)

        x = panel.xByModel{ii};
        y = panel.yByModel{ii};
        w = panel.wByModel{ii};
        m = panel.modelIndex(ii);

        if isempty(x) || numel(x) < minPoints
            continue
        end

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

        if smoothSig > 0
            H = gaussianSmooth2D(H, smoothSig);
        end

        lev = densityMassContourLevel(H, massLevel);

        if ~isfinite(lev) || lev <= 0
            continue
        end

        xc = 0.5 .* (xedges(1:end-1) + xedges(2:end));
        yc = 0.5 .* (yedges(1:end-1) + yedges(2:end));

        [~, h] = contour(ax, xc, yc, H.', [lev lev], ...
                         'LineColor', modelColors(m,:), ...
                         'LineWidth', lineWidth);

        set(h, 'HandleVisibility', 'off');

        hOut(ii) = h;
    end
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);

    ix = ix(keep);
    iy = iy(keep);
    w  = w(keep);

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

function lev = densityMassContourLevel(H, massLevel)

    z = H(:);
    z = z(isfinite(z) & z > 0);

    if isempty(z)
        lev = NaN;
        return
    end

    z = sort(z, 'descend');
    cs = cumsum(z) ./ sum(z);

    ii = find(cs >= massLevel, 1, 'first');

    if isempty(ii)
        lev = z(end);
    else
        lev = z(ii);
    end
end

function Hs = gaussianSmooth2D(H, sigma)

    if sigma <= 0
        Hs = H;
        return
    end

    rad = max(1, ceil(3*sigma));
    x = -rad:rad;

    g = exp(-(x.^2) ./ (2*sigma^2));
    g = g ./ sum(g);

    K = g' * g;

    Hs = conv2(H, K, 'same');
end

function stats = addWeightedLinearRegression(ax,x,y,w,xL,yL,confidence, ...
        lineColor,bandColor,bandAlpha,lineWidth,labelXFraction,fontSize)
%ADDWEIGHTEDLINEARREGRESSION Weighted OLS line, mean-response CI and R^2.
% Each model's biomass weights were normalized in collectPanelData, so the
% fit retains the figure's existing within-model biomass weighting while
% preventing mean biomass differences among models from dominating it.

    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);
    if numel(x) < 3 || range(x) <= eps
        return
    end

    % Probability weights and their Kish effective sample size.
    p = w ./ sum(w);
    nEff = (sum(w).^2) ./ sum(w.^2);
    if ~isfinite(nEff) || nEff <= 2
        return
    end

    X = [ones(size(x)), x];
    XtWX = X' * (X .* p);
    if rcond(XtWX) < 1e-12
        return
    end

    beta = XtWX \ (X' * (p .* y));
    yHat = X * beta;

    yBar = sum(p .* y);
    sse = sum(p .* (y-yHat).^2);
    sst = sum(p .* (y-yBar).^2);
    if sst > 0
        R2 = 1 - sse/sst;
    else
        R2 = NaN;
    end
    
    % Effective-sample-size-adjusted residual variance.
sigma2 = sse * nEff / (nEff-2);

% Covariance matrix of the fitted intercept and slope.
covBeta = (sigma2/nEff) .* inv(XtWX);

% Evaluate the regression across the observed x range.
xFitMin = max(xL(1),min(x));
xFitMax = min(xL(2),max(x));
xFit = linspace(xFitMin,xFitMax,250).';

XFit = [ones(size(xFit)),xFit];
yFit = XFit * beta;

% Standard error of the estimated mean regression.
seMean = sqrt(max(0,sum((XFit*covBeta).*XFit,2)));

% Confidence level is set by regressionConfidence in USER OPTIONS.
zCrit = sqrt(2) * erfinv(confidence);
ci = zCrit .* seMean;

lowerCI = yFit - ci;
upperCI = yFit + ci;

% Plot upper and lower confidence bounds.
% plot(ax,xFit,lowerCI,'--', ...
%     'Color',lineColor, ...
%     'LineWidth',1.2, ...
%     'HandleVisibility','off');
% 
% plot(ax,xFit,upperCI,'--', ...
%     'Color',lineColor, ...
%     'LineWidth',1.2, ...
%     'HandleVisibility','off');

% Plot the weighted regression over the confidence interval.
% plot(ax,xFit,yFit,'-', ...
%     'Color',lineColor, ...
%     'LineWidth',lineWidth, ...
%     'HandleVisibility','off');

    text(ax, 0.98, .48, sprintf('$R^2 = %.2f$',R2), ...
    'Units','normalized', ...
    'Interpreter','latex', ...
    'FontSize',fontSize, ...
    'FontWeight','bold', ...
    'Color',lineColor, ...
    'HorizontalAlignment','right', ...
    'VerticalAlignment','top', ...
    'Clipping','off', ...
    'HandleVisibility','off');

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

function addQuadPercW(ax, x, y, w)

    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);

    q = zeros(1,4);

    q(1) = 100*nansum(w(x <= 0 & y >  0)) / W;  % top-left
    q(2) = 100*nansum(w(x >  0 & y >  0)) / W;  % top-right
    q(3) = 100*nansum(w(x >  0 & y <= 0)) / W;  % bottom-right
    q(4) = 100*nansum(w(x <= 0 & y <= 0)) / W;  % bottom-left

    txtfmt = '\\textbf{%d\\%%}';
    m  = 0.015;
    fs = 18;

    text(ax, m, 1-m, sprintf(txtfmt, round(q(1))), ...
         'Units','normalized', ...
         'HorizontalAlignment','left', ...
         'VerticalAlignment','top', ...
         'FontSize',fs, ...
         'Interpreter','latex');

    text(ax, 1-m, 1-m, sprintf(txtfmt, round(q(2))), ...
         'Units','normalized', ...
         'HorizontalAlignment','right', ...
         'VerticalAlignment','top', ...
         'FontSize',fs, ...
         'Interpreter','latex');

    text(ax, 1-m, m, sprintf(txtfmt, round(q(3))), ...
         'Units','normalized', ...
         'HorizontalAlignment','right', ...
         'VerticalAlignment','bottom', ...
         'FontSize',fs, ...
         'Interpreter','latex');

    text(ax, m, m, sprintf(txtfmt, round(q(4))), ...
         'Units','normalized', ...
         'HorizontalAlignment','left', ...
         'VerticalAlignment','bottom', ...
         'FontSize',fs, ...
         'Interpreter','latex');
end

function placeSubplotLetters(axList, labels, offset, fs)

    if nargin < 3 || isempty(offset)
        offset = [-0.1 1];
    end

    if nargin < 4 || isempty(fs)
        fs = 18;
    end

    for k = 1:numel(axList)

        ax = axList(k);

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