%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Figure 1e - Scatter Plot %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Script Path
cd /Users/twrohr/Desktop/Matlab/Matlab/Grazing/CMIP6_Comp/NCC_Figs/Final_Figs/Revised/

% Variable Catalogue and Indecies 
for Variable_Catalogue = 1
%% Depth-Integrated Biomass Stocks
% 1  | Depth-Integrated Phytoplankton Biomass, 1–150 m (ΣP) | mmol m^-2
% 2  | Depth-Integrated Zooplankton Biomass, 1–150 m (ΣZ) | mmol m^-2
% 52 | Total Plankton Biomass (ΣBiomass = ΣP + ΣZ) | mmol m^-2

%% Bulk Carbon Fluxes
% 3  | Net Primary Production (NPP) | mmol m^-2 d^-1
% 5  | Export Production (EP) | mmol m^-2 d^-1
% 50 | Trophic Transfer (TT; mesozooplankton mortality to higher predation) | mmol m^-2 d^-1
% 4  | Grazing Losses (G) | mmol m^-2 d^-1
% 10 | Non-grazing Losses (NG) | mmol m^-2 d^-1
% 49 | Mesozooplankton Assimilation (GSPMeso; excludes microzooplankton) | mmol m^-2 d^-1
% 8  | Gross Secondary Production (GSP; all zooplankton, excludes assimilation inefficiencies) | mmol m^-2 d^-1
% 9  | Zooplankton Loss Rate (Z_m) | mmol m^-2 d^-1

%% Specific Rate Terms
% 6  | Phytoplankton-specific Division Rate (μ = NPP/ΣP) | d^-1
% 7  | Zooplankton-specific Grazing Rate (g = G/ΣZ) | d^-1
% 11 | Specific Grazing Loss Rate (l_g = G/ΣP) | d^-1
% 17 | Phytoplankton-specific Non-grazing Loss Rate (l_ng = NG/ΣP) | d^-1
% 16 | Phytoplankton-specific Loss Rate (l = l_g + l_ng) | d^-1

%% Community Composition
% 12 | Small Phytoplankton Composition (f_SP) | %
% 13 | Small Zooplankton Composition (f_SZ) | %
% 25 | Diatom Fraction | fraction
% 26 | Non-diatom Fraction | fraction

%% Physics
% 14 | Sea-surface Temperature (SST) | °C
% 15 | Mixed-layer Depth (MLD) | m

%% Carbon-Flux Transfer Efficiencies
% 18 | Export Efficiency (e-rat) | fraction
% 19 | Ecosystem Efficiency (eco-rat) | fraction
% 53 | Trophic Transfer Efficiency (TTE) | fraction

%% Phytoplankton Concentrations
% 23 | Zooplankton Biomass-weighted Phytoplankton Concentration ([Phy]_Zoo; concentration experienced by average zooplankton) | mmol m^-3
% 20 | Surface Phytoplankton Concentration ([Phyto]_Surf) | mmol m^-3
% 46 | Mean Phytoplankton Concentration within MLD ([Phyto]_MLD) | mmol m^-3

%% Vertical Profile
% 42 | Dilution Length Scale (DL) | m
% 21 | Deep Biomass Maximum Phytoplankton Concentration ([Phyto]_DBM) | mmol m^-3
% 22 | Deep Biomass Maximum Depth (Z_DBM) | m

%% Decomposed Loss-Rate Drivers: Grazing
% 27 | Direct Temperature Effect on Grazing Pressure (l_g,T) | d^-1
% 28 | Direct Vertical Redistribution Effect on Grazing Pressure (l_g,VD) | d^-1
% 29 | Direct Community Composition Effect on Grazing Pressure (l_g,CC) | d^-1
% 30 | Direct Phytoplankton Biomass Effect on Grazing Pressure (l_g,ΣP) | d^-1
% 31 | Direct Zooplankton Biomass Effect on Grazing Pressure (l_g,ΣZ) | d^-1

%% Decomposed Loss-Rate Drivers: Non-Grazing
% 32 | Direct Temperature Effect on Non-grazing Pressure (l_ng,T) | d^-1
% 33 | Direct Vertical Redistribution Effect on Non-grazing Pressure (l_ng,VD) | d^-1
% 34 | Direct Community Composition Effect on Non-grazing Pressure (l_ng,CC) | d^-1
% 35 | Direct Phytoplankton Biomass Effect on Non-grazing Pressure (l_ng,ΣP) | d^-1
% 36 | Direct Zooplankton Biomass Effect on Non-grazing Pressure (l_ng,ΣZ) | d^-1

%% Decomposed Loss-Rate Drivers: Total
% 37 | Direct Temperature Effect on Loss Rates (l_T) | d^-1
% 38 | Direct Vertical Redistribution Effect on Loss Rates (l_VD) | d^-1
% 39 | Direct Community Composition Effect on Loss Rates (l_CC) | d^-1
% 40 | Direct Phytoplankton Biomass Effect on Loss Rates (l_ΣP) | d^-1
% 41 | Direct Zooplankton Biomass Effect on Loss Rates (l_ΣZ) | d^-1
% 47 | Indirect Interaction Effects on Loss Rates (l_Int) | d^-1

%% Clearance Rates
% 44 | Clearance Rate (ε) | m^-3 mmol^-1 d^-1
% 45 | Direct Vertical Redistribution Effect on Clearance Rate (ε_VD) | m^-3 mmol^-1 d^-1
end

clear; clc;
close all;

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

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

% ---- Global look ----
set(groot, 'defaultAxesTickLabelInterpreter','latex', ...
           'defaultTextInterpreter','latex', ...
           'defaultLegendInterpreter','latex');

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

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

% ---- Variables ----
% Corrected to match labels:
% var 1 = Delta_% Sigma P / biomass
% var 6 = Delta_% division/growth rate
vX = 6;  %6
vY = 1;  %1

labelX = 'Change in phytoplankton growth rates; $\Delta_{\%} \mu\;(\%)$';
%labelX = '$\Delta_{\%} f_{sp}(\%)$';
labelY = 'Change in phytoplankton biomass; $\Delta_{\%} \Sigma \mathrm{P}\;(\%)$';

labelY = 'Change in Dilution Length Scale; $\Delta_{\%} D_L;(\%)$';
%labelX = '$\Delta_{\%} f_{sp}(\%)$';
labelX = 'Change in MLD; $\Delta_{\%} \; MLD (\%)$';


% ---- Plot settings ----
xmin = -80; xmax = 80; %125
ymin = -80; ymax = 80; %125
nbin = 250;            %250
 
% Percentile trimming; 0/100 = no trimming
outmin = 0;
outmax = 100;

% Same log-scaling idea as original plots.
% These are applied manually to the weighted bin density before converting
% each quadrant to RGB colour.
clL = 5e-6;
clU = 0.0014;

% Lower gamma = faint and mid-density structure becomes more visible.
density_gamma = .6;%.9;

% Slightly reduced boost so colours do not jump too quickly to the end-member.
contrast_boost = 1;%1.35;

% ---- Optional model contours on single ensemble panel ----
showModelContours = true;     % switch false to recover original no-contour look
contourMassLevel  = 0.85;     % contour enclosing densest 70% of model-weighted mass
contourNBins      = 200;
contourSmoothSig  = 3;      % smoothing in bin units; 0 = no smoothing
contourLW         = .7;
contourMinPoints  = 150;


% ---- Quadrant colour scheme ----
% Design:
%   white -> two shared cool greys -> three quadrant-specific colours
%
% This branches into colour earlier than the previous version, so more
% of the density field appears coloured rather than grey/black.


% Common low-density colours
white_rgb = [1.00, 1.00, 1.00];
cool1_rgb = [0.96, 0.97, 0.98];
cool2_rgb = [0.86, 0.88, 0.92];

% Common maximum-density colour
hot_rgb = [0,0,0];

% Bottom-left: blue branch
blue_1_rgb = [0.70, 0.83, 0.98];
blue_2_rgb = [0.43, 0.65, 0.90];
blue_3_rgb = [0.22, 0.45, 0.76];
blue_4_rgb = [0.08, 0.24, 0.55];

% Top-right: green branch
green_1_rgb = [0.82, 0.89, 0.57];
green_2_rgb = [0.60, 0.76, 0.35];
green_3_rgb = [0.38, 0.58, 0.19];
green_4_rgb = [0.18, 0.38, 0.08];

% Top-left: purple branch
purple_1_rgb = [0.86, 0.73, 0.94];
purple_2_rgb = [0.69, 0.50, 0.82];
purple_3_rgb = [0.51, 0.31, 0.68];
purple_4_rgb = [0.32, 0.12, 0.48];

% Bottom-right: orange branch
orange_1_rgb = [0.98, 0.76, 0.38];
orange_2_rgb = [0.91, 0.56, 0.16];
orange_3_rgb = [0.75, 0.37, 0.06];
orange_4_rgb = [0.52, 0.20, 0.02];

quadCols = struct();

% Keep the shared hot colour restricted to the highest densities
quadCols.knots = [ ...
    0.00 ... % white
    0.1 ... % pale grey
    0.2 ... % cool grey
    0.45 ... % quadrant colour 1
    0.7 ... % quadrant colour 2
    0.85 ... % quadrant colour 3
    0.95 ... % quadrant colour 4
    1];   % common hot endpoint

commonRamp = [ ...
    white_rgb; ...
    cool1_rgb; ...
    cool2_rgb];

quadCols.BL_ramp = [commonRamp; ...
    blue_1_rgb; blue_2_rgb; blue_3_rgb; blue_4_rgb; hot_rgb];

quadCols.TR_ramp = [commonRamp; ...
    green_1_rgb; green_2_rgb; green_3_rgb; green_4_rgb; hot_rgb];

quadCols.TL_ramp = [commonRamp; ...
    purple_1_rgb; purple_2_rgb; purple_3_rgb; purple_4_rgb; hot_rgb];

quadCols.BR_ramp = [commonRamp; ...
    orange_1_rgb; orange_2_rgb; orange_3_rgb; orange_4_rgb; hot_rgb];

% Labels for custom four-strip colour bars
quadBarLabels = {'','','',''};

% ---- Model contour 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],:) ...
];

modelColors = modelColors(1:11,:);

% ---- Optional weighted linear regression on ensemble panel ----
showRegression = true;
regLineColor   = [0 0 0];
regBandColor   = [0 0 0];
regLW          = 2.2;
regBandAlpha   = 0.006;
regTextFS      = 14;

%% =======================================================================
% Load model data once
% ========================================================================

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

Xmod = cell(1,11);
Ymod = cell(1,11);
Wmod = cell(1,11);

Xens = [];
Yens = [];
Wens = [];

nValidModels = 0;

for m = ms

    S = load(fullfile(dataDir, mods(mod_order(m)).name), ...
        'CoreVar_Inv_pc','Bio_Weights');

    x = squeeze(S.CoreVar_Inv_pc(vX,:,:));
    y = squeeze(S.CoreVar_Inv_pc(vY,:,:));
    w = S.Bio_Weights;

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

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

    if isempty(x)
        Xmod{m} = [];
        Ymod{m} = [];
        Wmod{m} = [];
        continue
    end

    % Model-normalization: each model total biomass weight = 1
    Wnorm = nansum(w);

    if ~isfinite(Wnorm) || Wnorm <= 0
        Xmod{m} = [];
        Ymod{m} = [];
        Wmod{m} = [];
        continue
    end

    w = w ./ Wnorm;

    % Optional trimming, model by model
    outsx = prctile(x,[outmin,outmax]);
    outsy = prctile(y,[outmin,outmax]);

    keep = x >= outsx(1) & x <= outsx(2) & ...
           y >= outsy(1) & y <= outsy(2) & ...
           isfinite(x) & isfinite(y) & isfinite(w) & w > 0;

    x = x(keep);
    y = y(keep);
    w = w(keep);

    Xmod{m} = x;
    Ymod{m} = y;
    Wmod{m} = w;

    Xens = [Xens; x];
    Yens = [Yens; y];
    Wens = [Wens; w];

    nValidModels = nValidModels + 1;
end

% Ensemble total weight approximately 1
if nValidModels > 0
    Wens = Wens ./ nValidModels;
end

%% =======================================================================
% Figure 1: single ensemble plot with model contours and marginal histograms
% ========================================================================

% Use fixed physical units so the main panel and marginal histograms align.
f1 = figure('Units','centimeters','Position',[2 2 28.2 26.3], 'Color','w');

% Axes positions: [left bottom width height], in cm
mainPos  = [2.7 2.6 23 23];
cbPos1 = [26 2.6 .8 23];  % [left bottom width height], cm
%topPos   = [2.7 23 20 1.85];
%rightPos = [23 2.6 1.85 20];

ax1     = axes('Parent',f1,'Units','centimeters','Position',mainPos);
%axTop   = axes('Parent',f1,'Units','centimeters','Position',topPos);
%axRight = axes('Parent',f1,'Units','centimeters','Position',rightPos);

% ---- Main quadrant density panel ----
plotQuadrantRGBDensity(ax1, Xens, Yens, Wens, nbin, [xmin xmax], [ymin ymax], ...
    clL, clU, density_gamma, contrast_boost, quadCols);

hold(ax1,'on');

% ---- Add individual model contours, no legend ----
if showModelContours
    plotModelContoursOnEnsemble(ax1, Xmod, Ymod, Wmod, modelColors, ...
        contourNBins, [xmin xmax], [ymin ymax], ...
        contourMassLevel, contourSmoothSig, contourLW, contourMinPoints);
end

% ---- Weighted linear regression with approximate 95% confidence band ----
if showRegression
    regStats = addWeightedLinearRegression(ax1, Xens, Yens, Wens, ...
        [xmin xmax], [ymin ymax], ...
        regLineColor, regBandColor, regLW, regBandAlpha, regTextFS);
end

% Redraw zero lines above contours
xline(ax1,0,'k','LineWidth',1.6);
yline(ax1,0,'k','LineWidth',1.6);

xlabel(ax1, labelX, 'FontSize', 25);
ylabel(ax1, labelY, 'FontSize', 25);

% Add quadrant percentages last so they sit above contours
addQuadPercW(ax1, Xens, Yens, Wens, 28);

% ---- Custom four-strip colour bar for RGB quadrant ramps ----
addQuadrantColorbars(f1, 'centimeters', cbPos1, quadCols, ...
    clL, clU, density_gamma, contrast_boost, quadBarLabels, 16);

% ---- Marginal histograms from same weighted 2D histogram ----
%plotMarginalHistograms(axTop, axRight, Xens, Yens, Wens, ...
%    nbin, [xmin xmax], [ymin ymax]);

%title(axTop, {'Ensemble'; 'all models, biomass-weighted, no TD/BU filter'}, ...
%    'Interpreter','latex','FontSize',20);

set(ax1, ...
    'FontSize', 28, ...
    'TickDir', 'in',  ...
    'TickLabelInterpreter','latex')
ax1.XAxis.FontSize = 28;
ax1.YAxis.FontSize = 28;



% =======================================================================
%% Helpers
% =======================================================================

function plotQuadrantRGBDensity(ax, x, y, w, nbin, xL, yL, clL, clU, density_gamma, contrast_boost, quadCols)
%PLOTQUADRANTRGBDENSITY
% One robust axes per panel. Computes weighted log-density and converts it
% to an RGB image with a different multi-stop ramp in each quadrant.
%
% Shared first half:
%   white -> cool greys -> charcoal
% Then each quadrant branches into its own colour family.

    cla(ax);

    if isempty(x) || isempty(y) || isempty(w)
        set(ax,'Visible','off');
        return
    end

    [Himg, xplot, yplot, Xc, Yc] = weightedHistogramImage(x, y, w, nbin, xL, yL);

    % Convert weighted density to log-scaled 0--1 intensity.
    logH = log10(Himg);
    t = (logH - log10(clL)) ./ (log10(clU) - log10(clL));
    t = max(min(t,1),0);

    % Stretch low-mid densities so the expanded early ramp is used.
    t = t .^ density_gamma;
    t = min(1, contrast_boost .* t);

    % Empty bins transparent
    alpha_img = double(isfinite(Himg));

    % RGB image starts white
    rgb = ones(size(Himg,1), size(Himg,2), 3);

    % Quadrants from bin centres
    qBL = Xc <= 0 & Yc <= 0;
    qTR = Xc >  0 & Yc >  0;
    qTL = Xc <= 0 & Yc >  0;
    qBR = Xc >  0 & Yc <= 0;

    rgb = applyDensityRamp(rgb, t, qBL, quadCols.BL_ramp, quadCols.knots);
    rgb = applyDensityRamp(rgb, t, qTR, quadCols.TR_ramp, quadCols.knots);
    rgb = applyDensityRamp(rgb, t, qTL, quadCols.TL_ramp, quadCols.knots);
    rgb = applyDensityRamp(rgb, t, qBR, quadCols.BR_ramp, quadCols.knots);

    h = imagesc(ax, xplot, yplot, rgb);
    set(h, 'AlphaData', alpha_img);

    set(ax, ...
        'YDir','normal', ...
        'Color',[1 1 1], ...
        'XLim',xL, ...
        'YLim',yL, ...
        'Box','on', ...
        'Layer','top', ...
        'FontSize',18, ...
        'LineWidth',1.1);

    axis(ax,'square');
    hold(ax,'on');

    xline(ax,0,'k','LineWidth',1.6);
    yline(ax,0,'k','LineWidth',1.6);
end


function [Himg, xplot, yplot, Xc, Yc] = weightedHistogramImage(x, y, w, nbin, xL, yL)
%WEIGHTEDHISTOGRAMIMAGE Weighted 2D histogram, image-ready orientation.

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

    H = accumarray([ix, iy], w, [nbin, nbin], @sum, 0);
    H(H <= 0) = NaN;

    % Transpose so rows correspond to y and columns to x
    Himg = H.';

    xplot = [xedges(1), xedges(end)];
    yplot = [yedges(1), yedges(end)];

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

function rgb = applyDensityRamp(rgb, t, qmask, ramp, knots)
%APPLYDENSITYRAMP
% Apply a multi-stop colour ramp to one quadrant.
%
% ramp  = nStop x 3 RGB values
% knots = 1 x nStop positions from 0 to 1
%
% Here the design is:
%   0.0 -> 0.5 : shared white / cool-grey / charcoal half
%   0.5 -> 1.0 : quadrant-specific colour half
%
% With knots:
%   [0, 1/6, 1/3, 1/2, 2/3, 5/6, 1]
%
% That gives equal colour resolution in the common and distinct halves.

    good = qmask & isfinite(t);

    if ~any(good(:))
        return
    end

    tt = t(good);

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

    rgb = max(min(rgb,1),0);
end



function plotModelContoursOnEnsemble(ax, Xmod, Ymod, Wmod, modelColors, ...
                                     nbin, xL, yL, massLevel, smoothSig, ...
                                     lineWidth, minPoints)
%PLOTMODELCONTOURSONENSEMBLE
% Adds one weighted density contour per model to the ensemble panel.
% No legend is created.

    hold(ax,'on');

    nModels = numel(Xmod);

    for m = 1:nModels

        x = Xmod{m};
        y = Ymod{m};
        w = Wmod{m};

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

        [Himg, ~, ~, Xc, Yc] = weightedHistogramImage(x, y, w, nbin, xL, yL);

        % Convert empty bins from NaN to 0 before smoothing/contouring
        H = Himg;
        H(~isfinite(H)) = 0;

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

        lev = densityMassContourLevel(H, massLevel);

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

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

        % No contour legend
        set(h, 'HandleVisibility', 'off');
    end
end


function lev = densityMassContourLevel(H, massLevel)
%DENSITYMASSCONTOURLEVEL
% Returns the density threshold enclosing the requested fraction of total
% weighted mass, starting from the densest bins.

    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)
%GAUSSIANSMOOTH2D Simple Gaussian smoothing without toolboxes.

    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 addQuadPercW(ax, x, y, w, fs)
%ADDQUADPERCW Add weighted quadrant percentages.

    if nargin < 5
        fs = 14;
    end

    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;  % TL
    q(2) = 100*nansum(w(x>0  & y>0))  / W;  % TR
    q(3) = 100*nansum(w(x>0  & y<=0)) / W;  % BR
    q(4) = 100*nansum(w(x<=0 & y<=0)) / W;  % BL

    txtfmt = '\\textbf{%d\\%%}';
    m = 0.005; m1 = .395; m2 = .395;

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

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

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

    text(ax, m1, 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


function plotMarginalHistograms(axTop, axRight, x, y, w, nbin, xL, yL)
%PLOTMARGINALHISTOGRAMS
% Adds top and right marginal weighted histograms using the same 2D bins as
% the main quadrant density plot.
%
% Top histogram:
%   integrates across y / rows, showing weighted distribution along x.
%
% Right histogram:
%   integrates across x / columns, showing weighted distribution along y.

    [Himg, ~, ~, Xc, Yc] = weightedHistogramImage(x, y, w, nbin, xL, yL);

    H = Himg;
    H(~isfinite(H)) = 0;

    xcent = Xc(1,:);
    ycent = Yc(:,1);

    hx = nansum(H,1);      % x marginal: integrate over y / rows
    hy = nansum(H,2);      % y marginal: integrate over x / columns

    histCol = [0.25 0.25 0.25];

    % ---- Top marginal histogram ----
    cla(axTop);
    bar(axTop, xcent, hx, 1.0, ...
        'FaceColor',histCol, ...
        'EdgeColor','none', ...
        'FaceAlpha',0.75);

    hold(axTop,'on');
    xline(axTop,0,'k','LineWidth',1.2);

    xlim(axTop,xL);

    if max(hx) > 0
        ylim(axTop,[0 max(hx).*1.12]);
    end

    set(axTop, ...
        'XTickLabel',[], ...
        'TickLabelInterpreter','latex', ...
        'FontSize',18, ...
        'Box','on', ...
        'Layer','top', ...
        'LineWidth',1.0);

    ylabel(axTop,'Weighted freq.', ...
        'Interpreter','latex', ...
        'FontSize',11);

    % ---- Right marginal histogram ----
    cla(axRight);
    barh(axRight, ycent, hy, 1.0, ...
        'FaceColor',histCol, ...
        'EdgeColor','none', ...
        'FaceAlpha',0.75);

    hold(axRight,'on');
    yline(axRight,0,'k','LineWidth',1.2);

    ylim(axRight,yL);

    if max(hy) > 0
        xlim(axRight,[0 max(hy).*1.12]);
    end

    set(axRight, ...
        'YTickLabel',[], ...
        'TickLabelInterpreter','latex', ...
        'FontSize',18, ...
        'Box','on', ...
        'Layer','top', ...
        'LineWidth',1.0);

    xlabel(axRight,'Weighted freq.', ...
        'Interpreter','latex', ...
        'FontSize',11);
end


function stats = addWeightedLinearRegression(ax, x, y, w, xL, yL, ...
                                             lineColor, bandColor, ...
                                             lineWidth, bandAlpha, fs)
%ADDWEIGHTEDLINEARREGRESSION
% Adds weighted linear regression line, visible approximate 95% confidence
% lines, optional confidence band, and weighted R2 annotation.
%
% Fit is restricted to 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
        warning('Not enough valid points for weighted regression.');
        return
    end

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

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

    if neff <= p + 1
        warning('Effective sample size too small for regression uncertainty.');
        return
    end

    % Weighted least squares: y = b0 + b1*x
    X = [ones(size(x)), x];

    XtWX = X' * (X .* w);
    XtWy = X' * (y .* w);

    beta = XtWX \ XtWy;

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

    % Weighted R2
    ybar = nansum(w .* y);
    SSE  = nansum(w .* resid.^2);
    SST  = nansum(w .* (y - ybar).^2);
    R2   = 1 - SSE ./ SST;

    % Approximate residual variance.
    % This treats the normalized weights as probability weights and uses
    % effective sample size for uncertainty.
    sigma2 = SSE .* neff ./ max(neff - p, 1);

    % Coefficient covariance
    covB = sigma2 .* inv(XtWX);

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

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

    slopeCI95 = tcrit .* slopeSE;

    % Prediction grid
    xg = linspace(xL(1), xL(2), 300).';
    Xg = [ones(size(xg)), xg];

    yg = Xg * beta;

    % 95% confidence interval for mean fitted line
    varMean = sum((Xg * covB) .* Xg, 2);
    ci = tcrit .* sqrt(max(varMean, 0));

    ylo = yg - ci;
    yhi = yg + ci;

    hold(ax,'on');

    % Optional confidence band. Do NOT push this to the bottom or it can
    % disappear behind the density image.
    if bandAlpha > 0
        fill(ax, ...
            [xg; flipud(xg)], ...
            [ylo; flipud(yhi)], ...
            bandColor, ...
            'EdgeColor','none', ...
            'FaceAlpha',bandAlpha, ...
            'HandleVisibility','off');
    end

    % Explicit confidence boundary lines
    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');

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

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


function cbAx = addQuadrantColorbars(parentFig, units, pos, quadCols, ...
                                     clL, clU, density_gamma, ...
                                     contrast_boost, quadLabels, fs)
%ADDQUADRANTCOLORBARS
% Custom RGB colour bars for quadrant-specific density ramps.
%
% This is needed because the main panels are RGB images, so MATLAB's
% standard colorbar does not know about the quadrant-specific colour maps.
%
% Four adjacent vertical strips are shown:
%   BL = x <= 0, y <= 0
%   TR = x >  0, y >  0
%   TL = x <= 0, y >  0
%   BR = x >  0, y <= 0
%
% The vertical axis is biomass-weighted frequency, matching the same
% log scaling, gamma, and contrast boost used in the main panels.

    if nargin < 10 || isempty(fs)
        fs = 10;
    end

    n = 256;

    logVals = linspace(log10(clL), log10(clU), n).';
    Hvals   = 10.^logVals;

    % Same transformation used in plotQuadrantRGBDensity
    t = (log10(Hvals) - log10(clL)) ./ (log10(clU) - log10(clL));
    t = max(min(t,1),0);

    t = t .^ density_gamma;
    t = min(1, contrast_boost .* t);

    ramps = { ...
        quadCols.BL_ramp, ...
        quadCols.TR_ramp, ...
        quadCols.TL_ramp, ...
        quadCols.BR_ramp};

    rgbBar = ones(n,4,3);

    for q = 1:4
        ramp = ramps{q};

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

    rgbBar = max(min(rgbBar,1),0);

    cbAx = axes('Parent', parentFig, ...
        'Units', units, ...
        'Position', pos);

    imagesc(cbAx, 1:4, logVals, rgbBar);
    set(cbAx, ...
        'YDir','normal', ...
        'YAxisLocation','right', ...
        'XTick',[], ...
        'XTickLabel',[], ...
        'TickLabelInterpreter','latex', ...
        'FontSize',fs, ...
        'Box','on', ...
        'LineWidth',1.0, ...
        'Layer','top');

    xlim(cbAx,[0.5 4.5]);
    ylim(cbAx,[log10(clL) log10(clU)]);

    % Nice log tick labels
    tickVals = [1e-5 1e-4 1e-3 1e-2];
    tickVals = tickVals(tickVals >= clL & tickVals <= clU);

    tickLabs = cell(size(tickVals));
    for i = 1:numel(tickVals)
        tickLabs{i} = sprintf('$10^{%d}$', round(log10(tickVals(i))));
    end

    set(cbAx, ...
        'YTick',log10(tickVals), ...
        'YTickLabel',tickLabs);

    ylabel(cbAx, '', ...
        'Interpreter','latex', ...
        'FontSize',fs+1);

    xlabel(cbAx, '', ...
        'Interpreter','latex', ...
        'FontSize',fs);

    title(cbAx, '', ...
        'Interpreter','latex', ...
        'FontSize',fs+1);

    hold(cbAx,'on');

    % Separators between the four adjacent strips
    for xx = .5:1:3.5
        plot(cbAx, [xx xx], [log10(clL) log10(clU)], '-', ...
            'Color',[0.25 0.25 0.25], ...
            'LineWidth',0, ...
            'HandleVisibility','off');
    end
end