%%%%%%%%%%%%%%%%%%%%%%%%%
%% Fig. 4b-g — Ensemble Maps of loss rate Driver Terms
%%%%%%%%%%%%%%%%%%%%%%%%%
cd /Users/twrohr/Desktop/Matlab/Matlab/Grazing/CMIP6_Comp/NCC_Figs/Final_Figs/Revised
addpath(genpath('/Users/twrohr/Desktop/Matlab/Matlab/Utilities'))
% IMPORTANT: td_frequency from Fig2_AllMaps_exp must exist
% should these be 

for startup =1  
% Load regridded lat lon
cd /Users/twrohr/Desktop/data.nosync/SSP585
         load('Longhurst.mat');
          Lon =  LonGrid; Lat = LatGrid;

% Constants
Ra = 6.371e6;                  % Earth radius [m]
deg2rad = pi / 180;
dphi = deg2rad;               % 1 degree in radians
dlambda = deg2rad;

% Latitude vector (center of cells)
lat = linspace(-89.5, 89.5, 180);  % 180 latitudes from -89.5 to +89.5
lat_rad = lat * deg2rad;

% Area per 1x1 degree cell as a function of latitude
area_per_lat = Ra^2 * dphi * dlambda * cos(lat_rad);  % [m^2]

% Repeat across 360 longitudes to make full 180x360 grid
area_map = repmat(area_per_lat', 1, 360);  % [180x360]

% Convert lon for warppign --- Shift longitudes to [0, 360]
        lonvec = linspace(-179.5, 179.5, 360);
        latvec = linspace(-89.5, 89.5, 180);
        lonvec_wrapped = lonvec;
        lonvec_wrapped(lonvec_wrapped < 0) = lonvec_wrapped(lonvec_wrapped < 0) + 360;
        [lonvec_sorted, lon_order] = sort(lonvec_wrapped);

% Colormaps
    % Get full bam colormap (e.g., 256 colors)
    n_colors = 256;

    % --- Define BU Map  ---
    blue_rgb   = [0.35, 0.55, 0.85];  % muted slate blue
    white_rgb  = [1.0,  1.0,  1.0];   % white
    green_rgb  = [0.30, 0.65, 0.35];  % classy muted green
        % --- Interpolation ---
        n_half = floor(n_colors/2);
        blue_to_white = [linspace(blue_rgb(1), white_rgb(1), n_half)', ...
                         linspace(blue_rgb(2), white_rgb(2), n_half)', ...
                         linspace(blue_rgb(3), white_rgb(3), n_half)'];
        white_to_green = [linspace(white_rgb(1), green_rgb(1), n_colors - n_half)', ...
                          linspace(white_rgb(2), green_rgb(2), n_colors - n_half)', ...
                          linspace(white_rgb(3), green_rgb(3), n_colors - n_half)'];
        % --- Combine into full colormap ---
        cmapBU = [blue_to_white; white_to_green];
    
    % --- Define TD map  ---
    red_rgb     = [0.78, 0.25, 0.25];  % soft classy red
    mustard_rgb = [0.90, 0.7, 0.25];   % strong yellow
        % --- Interpolation ---
        red_to_white = [linspace(red_rgb(1), white_rgb(1), n_half)', ...
                        linspace(red_rgb(2), white_rgb(2), n_half)', ...
                        linspace(red_rgb(3), white_rgb(3), n_half)'];
        white_to_mustard = [linspace(white_rgb(1), mustard_rgb(1), n_colors - n_half)', ...
                            linspace(white_rgb(2), mustard_rgb(2), n_colors - n_half)', ...
                            linspace(white_rgb(3), mustard_rgb(3), n_colors - n_half)'];
        % --- Combine into full colormap ---
        cmapTD = [red_to_white; white_to_mustard];

    %%% Discrete
    green_start = [0.40, 0.60, 0.40];
    green_end   = [0.85, 0.95, 0.85];
    greens = [linspace(green_start(1), green_end(1), 6)', ...
              linspace(green_start(2), green_end(2), 6)', ...
              linspace(green_start(3), green_end(3), 6)'];
    red_start = [0.90, 0.70, 0.70];
    red_end   = [0.70, 0.00, 0.00];
    reds = [linspace(red_start(1), red_end(1), 6)', ...
            linspace(red_start(2), red_end(2), 6)', ...
            linspace(red_start(3), red_end(3), 6)'];
    custom_cmap = [greens; reds];

    % Model Colors
    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,12];
    cmap1 = brewermap(13,'Oranges');
    cmap2 = brewermap(15,'Greens');
    cmap3 = brewermap(11,'Purples');
    cmap_mod = [cmap2([5,7,11,13,15],:); cmap1([5,7,11,13],:); cmap3([7,11],:); 0 0 0];

    CoreVar_Names = {'Phytoplankton Biomass',     'Zooplankton Biomass',         'NPP',                        'Grazing Losses',            'EP',                 ...
                     'Spec. Division Rate',       'Spec. Grazing Rate',          'GSP',                        'Z Loss rate',               'Non-Grazing Losses', ...
                     'Spec. Grazing Loss Rate',   'Fraction Small P',            'Fraction Small Z',           'SST' ,                      'MLD' ,               ...
                     'Spec. Loss Rate',           'Spec. Non-grazing Loss Rate'  'Export Efficeincy'           'Ecosystem Efficeincy'       'Surface [Phyto]' ,  ...
                     'DBM [Phyto]',               'DBM Depth',                   'Dilution Index',             'Diaz Frac',                 'Diatom Frac',           'Non- Diatom Frac'...
                     'Temperature (g)',           'Vertical Dist (g) ',          'Community Composition (g)',  'Phyto. Biomass (g)',        'Zooplankton Biomass (g)',...
                     'Temperature (n)',           'Vertical Dist (n) ',          'Community Composition (n)',  'Phyto. Biomass (n)',        'Zooplankton Biomass (n)',...
                     'Temperature (T)',           'Vertical Dist (T) ',          'Community Composition (T)',  'Phyto. Biomass (T)',        'Zooplankton Biomass (T)',...
                     'Dilution Index',};

     decomp_names = {{'$\Delta l_{T}$'; 'First-order Effect of' ; 'Temperture'}, {'$\Delta l_{VD}$'; 'First-order Effect of' ; 'Vertical Distribution'},...
                     {'$\Delta l_{CC}$';'First-order Effect of' ; 'Community Composition'}, {'$\Delta l_{\Sigma P}$';'First-order Effect of' ; 'Phytoplankton Biomass'},...
                     {'$\Delta l_{\Sigma Z}$';'First-order Effect of';'Zooplankton Biomass'}, {'$\Delta l_{I}$';'First-order Effect of';'Interaction Effects'}};
     decomp_symbols = {{'$\Delta l_{T}$ (d$^{-1}$)'}, {'$\Delta l_{VD}$ (d$^{-1}$)'},...
                     {'$\Delta l_{CC}$ (d$^{-1}$)'}, {'$\Delta l_{\Sigma P}$ (d$^{-1}$)'},...
                     {'$\Delta l_{\Sigma Z}$ (d$^{-1}$)'},{'$\Delta l_{Int}$ (d$^{-1}$)'}};
    
     decomp_names2 =  { '$\frac{\Delta \; l_{T}}{\Delta\;  \mathrm{SST}}$' , '$\frac{\Delta \; l_{VD} }{\Delta \; \mathrm{D}_L}$' ,...
                    '$\frac{\Delta \; l_{CC} }{\Delta\;  f_{\mathrm{SP}}}$' , '$\frac{\Delta \; l_{\Sigma P}}{\Delta \; \Sigma\,{\mathrm{P}}}$' ,...
                     '$\frac{\Delta \; l_{\Sigma Z} }{\Delta\;  \Sigma \mathrm{P} }$' };

     decomp_units = {'(d$^{-1}$ $^{\circ}$C warmer $^{-1}$)', '(d$^{-1}$ m $^{-1}$)', '(d$^{-1}$ $\%$ inc. small phy. $^{-1}$)', '(m$^{2}$ d$^{-1}$ mmol$^{-1}$)', '(m$^{2}$ d$^{-1}$ mmol$^{-1}$)'};

     Driver_label = {'$\Delta$ SST ($^\circ$C)',' $\Delta$ Dilution Length Scale ($D_{L}$; m)', ' $\Delta$ Fraction Sm. Phytoplankton ($f_{sp}$; $\%$)', '$\Delta$ $\Sigma$ Phytoplankton (mmolC m$^{-2}$)',...
                     '$\Delta$ $\Sigma$ Zooplankton (mmolC m$^{-2}$)','dummy'};

     % decomp_names =decomp_names([4,5,1,3,2]);
     % decomp_names2 =decomp_names2([4,5,1,3,2]);
     % Driver_label =Driver_label([4,5,1,3,2]);
     % decomp_units =decomp_units([4,5,1,3,2]);
     % decomp_symbols =decomp_symbols([4,5,1,3,2]); 
     CoreVar_Names{37} = '-';

     wedge_colors = [
    0.90 0.55 0.25;   % SST     → soft orange
    0.35 0.60 0.85;   % MLD     → clearer blue
    0.25 0.75 0.65;   % Diatoms → more greenish teal
    0.60 0.75 0.45;   % Phyto   → sage green
    0.85 0.55 0.75    % Zoo     → soft magenta
];
lighter_colors = wedge_colors + 0.6 .* (1 - wedge_colors);  % shift toward white, preserve hue

end

%% === Loss Rate: six maps in a 2 x 3 layout =============================

% =========================================
clc; close all
clearvars plot

% ---------------- CONFIG ----------------
% Panel order (left-to-right, then top-to-bottom):
% Phytoplankton biomass, Zooplankton biomass, Temperature,
% Community composition, Vertical distribution, Interaction effects.
vord        = [4,5,1,3,2,6];
VARIDS      = [37,38,39,40,41,47];   % 37,38,39,40,41,47];
var_ids     = VARIDS(vord);        
vars2plot   = 1:6;
panel_titles = {'Phytoplankton Biomass', 'Zooplankton Biomass', ...
                'Temperature', 'Community Composition', ...
                'Vertical Distribution', 'Interaction Effects'};

% Dynamic latitude window — set once
minLat     = -65;        % <<< adjust as desired
maxLat     =  65;        % <<< adjust as desired

% ---------------- DATA ----------------
basedir    = '/Users/twrohr/Desktop/data.nosync/SSP585';
load(fullfile(basedir,'All_Stats_Raw.mat'));   % expects: InterMod_mean_maps, Lat_t, Lon_t, latvec, zon_stats, cmap_mod, etc.
%load(fullfile(basedir,'All_Stats_Norm.mat'));  % For partials normalized by change in Delta l

% ---------------- 2 x 3 MAP LAYOUT ----------------
nPanelRows = 2;
nPanelCols = 3;
left_m     = 0.025;
right_m    = 0.025;
top_m      = 0.055;
bot_m      = 0.105;       % space for the shared colorbar
col_gap    = 0.025;
row_gap    = 0.055;

panel_w = (1 - left_m - right_m - (nPanelCols-1)*col_gap) / nPanelCols;
panel_h = (1 - top_m - bot_m - (nPanelRows-1)*row_gap) / nPanelRows;

fig = figure('Units','normalized','Position',[0.02 0.08 0.96 0.82], ...
             'Color','w');
ax_maps = gobjects(numel(vars2plot),1);

% =========================================================
% ---------------- LOOP OVER SIX MAP PANELS ----------------
% =========================================================
for r = vars2plot
    pidx   = var_ids(r);

    panel_row = ceil(r / nPanelCols);
    panel_col = mod(r-1, nPanelCols) + 1;
    x0 = left_m + (panel_col-1) * (panel_w + col_gap);
    y0 = 1 - top_m - panel_row*panel_h - (panel_row-1)*row_gap;

    pos_map = [x0, y0, panel_w, panel_h];
    axm = axes('Units','normalized','Position',pos_map, ...
               'PositionConstraint','innerposition');
    ax_maps(r) = axm;
    axes(axm);
    axesm('robinson','MapLonLimit',[0 360], ...
          'MapLatLimit',[minLat maxLat],'Origin',[0 180 0]);
    axis off; tightmap

    % Inter-model mean for this variable (TOTALS: vidx+10), mask to lat window
    plotvar = squeeze(InterMod_mean_maps(pidx,:,:)); % 
    lat_rows = Lat_t(:,1);
    mask_lat = (lat_rows < minLat) | (lat_rows > maxLat);
    plotvar(mask_lat,:) = NaN;

    Ra = georefcells([min(Lat_t(:,1)) max(Lat_t(:,1))], ...
                     [min(Lon_t(1,:)) max(Lon_t(1,:))], size(plotvar));
    geoshow(plotvar, Ra, 'DisplayType','texturemap');
    load coastlines; patchm(coastlat,coastlon,[.7 .7 .7]);

    % colors (reversed balance)
    % Use fixed or auto-centered limits as before
    climv = [ -.062, .062];
    %climv = [ -150, 150]; % For Normalized partials (divided by total loss change)
    %climv = [ -50, 50]; % For  SP fraction
    colormap(axm, flipud(cmocean('balance')));
    clim(axm,climv);

    title(axm, panel_titles{r}, 'FontSize',18, 'FontWeight','normal', ...
          'Interpreter','latex');


end % rows

% One shared horizontal colorbar for all six panels.
last_map_pos = get(ax_maps(end),'Position');
cb = colorbar(ax_maps(end),'southoutside');
cb.Units = 'normalized';
cb.Position = [0.33 0.045 0.34 0.024];
cb.TickLabelInterpreter = 'latex';
cb.FontSize = 13;
cb.Label.String = '$\Delta l$ (d$^{-1}$)';
cb.Label.Interpreter = 'latex';
cb.Label.FontSize = 15;
set(ax_maps(end),'Position',last_map_pos);


%% =========================================================
% Correlation between driver maps and td_frequency
% =========================================================
% Computes Pearson correlation between the spatial distribution of each
% plotted driver map and td_frequency, restricted to 65S--65N.
%
% Reports:
%   1) full correlation over all valid cells
%   2) correlation where driver variable is increasing, D > 0
%   3) correlation where driver variable is decreasing, D < 0
%
% Requires td_frequency to already be in workspace, as noted at top of script.

fprintf('\n============================================================\n')
fprintf('Spatial Pearson correlations with td_frequency, 65S--65N\n')
fprintf('Full, driver-increasing cells, and driver-decreasing cells\n')
fprintf('============================================================\n')

if ~exist('td_frequency','var')
    error('td_frequency is not in the workspace. Run the script that creates td_frequency first.')
end

% Ensure td_frequency is 2D
E = squeeze(td_frequency);

% Latitude mask: 65S to 65N
lat_mask_65 = Lat_t >= -65 & Lat_t <= 65;

% If needed, try to handle transposed td_frequency
if ~isequal(size(E), size(Lat_t))
    if isequal(size(E'), size(Lat_t))
        E = E';
        warning('td_frequency was transposed to match Lat_t.')
    else
        error('td_frequency size does not match Lat_t. Check dimensions of td_frequency, Lat_t, and InterMod_mean_maps.')
    end
end

% Area weights, if available. If not, use unweighted Pearson.
use_area_weighting = exist('area_map','var') && isequal(size(area_map), size(Lat_t));

if use_area_weighting
    corr_type = 'area-weighted Pearson r';
else
    corr_type = 'unweighted Pearson r';
end

fprintf('Using: %s\n\n', corr_type)

for rr = vars2plot

    pidx = var_ids(rr);

    % This is the same map variable plotted in the first-column map panel
    D = squeeze(InterMod_mean_maps(pidx,:,:));

    % Basic dimension check
    if ~isequal(size(D), size(E))
        if isequal(size(D'), size(E))
            D = D';
            warning('Driver map pidx=%d was transposed to match td_frequency.', pidx)
        else
            error('Driver map pidx=%d size does not match td_frequency.', pidx)
        end
    end

    valid_base = lat_mask_65 & isfinite(D) & isfinite(E);

    if use_area_weighting
        W = area_map;
        valid_base = valid_base & isfinite(W) & W > 0;

        valid_all = valid_base;
        valid_inc = valid_base & D > 0;
        valid_dec = valid_base & D < 0;

        r_all = local_weighted_corr_simple(D(valid_all), E(valid_all), W(valid_all));
        r_inc = local_weighted_corr_simple(D(valid_inc), E(valid_inc), W(valid_inc));
        r_dec = local_weighted_corr_simple(D(valid_dec), E(valid_dec), W(valid_dec));

    else
        valid_all = valid_base;
        valid_inc = valid_base & D > 0;
        valid_dec = valid_base & D < 0;

        r_all = corr(D(valid_all), E(valid_all), 'Type','Pearson', 'Rows','complete');
        r_inc = corr(D(valid_inc), E(valid_inc), 'Type','Pearson', 'Rows','complete');
        r_dec = corr(D(valid_dec), E(valid_dec), 'Type','Pearson', 'Rows','complete');
    end

    % Make readable label
    if exist('decomp_symbols','var') && numel(decomp_symbols) >= vord(rr)
        label = decomp_symbols{vord(rr)};
        if iscell(label)
            label = strjoin(label(:)', ' ');
        end
    elseif exist('decomp_names','var') && numel(decomp_names) >= vord(rr)
        label = decomp_names{vord(rr)};
        if iscell(label)
            label = strjoin(label(:)', ' ');
        end
    else
        label = sprintf('pidx %d', pidx);
    end

    fprintf('%-45s | pidx = %3d\n', label, pidx)
    fprintf('   all cells      : n = %6d | r = %+0.3f\n', nnz(valid_all), r_all)
    fprintf('   driver > 0     : n = %6d | r = %+0.3f\n', nnz(valid_inc), r_inc)
    fprintf('   driver < 0     : n = %6d | r = %+0.3f\n', nnz(valid_dec), r_dec)
    fprintf('\n')
end

fprintf('============================================================\n\n')

% =========================
% ----- Helpers ----------
% =========================
function out = ternary(cond, a, b)
if cond, out = a; else, out = b; end
end

function counts = histcounts_weighted(x, edges, w)
% Weighted histogram per bin edge
    x = x(:); w = w(:);
    idx   = discretize(x, edges);
    ok    = ~isnan(idx) & isfinite(w);
    counts_col = accumarray(idx(ok), w(ok), [numel(edges)-1, 1], @sum, 0);
    counts = counts_col.';   % row
end

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

function out = smartfmt(v)
% Always scientific notation, 2 significant figures, 2 lines
    if ~isfinite(v) || v == 0
        mant = 0;  exp10 = 0;
    else
        exp10 = floor(log10(abs(v)));
        mant  = v / 10^exp10;
        mant  = round(mant, 2, 'significant');
        if abs(mant) >= 10
            mant  = mant / 10;
            exp10 = exp10 + 1;
        end
    end
    mant_str = sprintf('%+.1f', mant);
    if mant == 0 && exp10 ==0
        line1 = ''; line2 = '';
    else
        line1 = ['$' mant_str '$'];
        line2 = ['e ' num2str(exp10)];
    end
    out = {line1; line2};
end

function c = signcolor(v)
% Blue if +, red if -, gray if ~0/NaN
    if ~isfinite(v) || v==0
        c = [0.35 0.35 0.35];
    elseif v > 0
        c = [0.10 0.45 0.90];
    else
        c = [0.80 0.20 0.20];
    end
end

function draw_half_box(ax, c, side, med, q25, q75, maxW, col)
    if ~all(isfinite([med,q25,q75])) || q75<=q25, return; end
    w  = 0.55*maxW;
    x0 = c + strcmp(side,'left')*(-w) + strcmp(side,'right')*(0);
    r  = rectangle(ax,'Position',[x0, q25, w, (q75-q25)], ...
        'FaceColor',[0 0 0 0.12],'EdgeColor',col,'LineWidth',0.8);
    m  = plot(ax, [x0, x0+w], [med, med], 'Color', col, 'LineWidth', 1.6);
    uistack([r m],'top');
end

function f = weighted_frac(x, w, cond)
    valid = isfinite(x) & isfinite(w) & (w > 0);
    x    = x(valid);
    w    = w(valid);
    cond = cond(valid);

    den = sum(w);
    if den == 0
        f = NaN;
    else
        f = sum(w(cond)) / den;
    end
end

function r = local_weighted_corr_simple(x, y, w)
% Area-weighted Pearson correlation coefficient.

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

    x = x(ok);
    y = y(ok);
    w = w(ok);

    if numel(x) < 3 || sum(w) <= 0
        r = NaN;
        return
    end

    w = w(:) ./ sum(w);
    x = x(:);
    y = y(:);

    mx = sum(w .* x);
    my = sum(w .* y);

    dx = x - mx;
    dy = y - my;

    vx = sum(w .* dx.^2);
    vy = sum(w .* dy.^2);

    if vx <= 0 || vy <= 0
        r = NaN;
    else
        r = sum(w .* dx .* dy) ./ sqrt(vx .* vy);
    end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%% End — Combined Final Script %%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
