plot_hierarchical_results.m

function plot_hierarchical_results
% Plots for hierarchical IMR comparison.
% Expects 'hierarchical_results.mat' with: Models, post, logZ_vec, modelNames, betaGrid.

clear; clc;

modelNames = {'Newt','NH','NHKV','qNH','LM','qKV','SLS'};

matFile = 'hierarchical_results.mat';
assert(isfile(matFile), 'File not found: %s', matFile);
S = load(matFile);

% ---- sanity checks ----
req = {'Models','post','logZ_vec','modelNames','betaGrid'};
for k = 1:numel(req)
    assert(isfield(S, req{k}), 'Missing "%s" in %s', req{k}, matFile);
end

Models     = S.Models;
post       = S.post(:);
logZ_vec   = S.logZ_vec(:);
% modelNames = S.modelNames(:);
betaGrid   = S.betaGrid(:);

outDir = fileparts(which(matFile));
if isempty(outDir), outDir = pwd; end

% ---- style pools (markers + colors) ----
markers = {'o','s','d','^','v','>','<','p','h','x','+'};  % wraps if needed
markerSizeBig = 12;  % big markers
lineThick     = 3;   % thick lines

% Custom color set avoiding yellow
colors = [
    0.00 0.45 0.74;  % blue
    0.85 0.33 0.10;  % reddish-brown
    0.47 0.67 0.19;  % green
    0.49 0.18 0.56;  % purple
    0.30 0.75 0.93;  % cyan
    0.64 0.08 0.18;  % dark red
    0.00 0.00 0.55;  % navy
    0.35 0.35 0.35;  % dark gray
];
if size(colors,1) < numel(modelNames)
    colors = repmat(colors, ceil(numel(modelNames)/size(colors,1)), 1);
end

% ---------- 1) Model posterior bar chart ----------
fig1 = figure('Color','w');
bar(post, 'FaceColor',[0.3 0.6 0.9], 'EdgeColor','none');
xticks(1:numel(modelNames));
xticklabels(modelNames);
ylabel('$P(M\mid D)$','Interpreter','latex');
ylim([0 1]);
set(gca,'XGrid','off','YGrid','off');
formatPlot();
savepng(fig1, fullfile(outDir, 'fig_model_posteriors.png'));

% ---------- 2) Log-evidence (logZ) comparison ----------
fig2 = figure('Color','w');
bar(logZ_vec, 'FaceColor',[0.5 0.5 0.5], 'EdgeColor','none');
xticks(1:numel(modelNames));
xticklabels(modelNames);
ylabel('$\log_{10} P(D\mid M)$','Interpreter','latex');
set(gca,'XGrid','off','YGrid','off');
formatPlot();
savepng(fig2, fullfile(outDir, 'fig_logZ_comparison.png'));

% ---------- 3) β posterior for the winning model ----------
% ---------- 4) β posterior comparison across models (viridis + proper hold) ----------
fig4 = figure('Color','w');
ax4  = axes('Parent',fig4);

% Viridis palette (fallback to parula if viridis.m isn't on path)
if exist('viridis','file')==2
    C = viridis(numel(modelNames) + 2);
else
    C = parula(numel(modelNames) + 2);
end
C = C(2:end-1,:);                       % trim extremes for nicer mid-tones
set(ax4, 'ColorOrder', C);              % set palette
hold(ax4, 'on');                        % <-- keep all subsequent plots

peakIdx = [];                           % rightmost peak across models
h = gobjects(0);                        % for legend handles

for m = 1:numel(modelNames)
    bp = Models(m).post_beta(:);
    if numel(bp)==numel(betaGrid) && all(isfinite(bp))
        % advance ColorOrder by model index (optional; MATLAB cycles automatically)
        set(ax4,'ColorOrderIndex', m);
        h(end+1) = plot(betaGrid, bp, '-', ...
            'LineWidth', lineThick, ...
            'DisplayName', modelNames{m}); %#ok<AGROW>
        [~,pk] = max(bp);
        peakIdx(end+1) = pk; %#ok<AGROW>
    end
end

% x-limit: show up to a bit past the rightmost peak among all models
if ~isempty(peakIdx)
    right = min(betaGrid(end), betaGrid(max(peakIdx)) + 0.25);
    xlim([betaGrid(1) right]);
end

xlabel('$\beta$','Interpreter','latex');
ylabel('$P(\beta\mid D,M)$','Interpreter','latex');
legend(h, 'Location','bestoutside','Interpreter','latex');
set(ax4,'XGrid','off','YGrid','off');
formatPlot();
savepng(fig4, fullfile(outDir, 'fig_beta_posteriors_all_models.png'));

% ---------- 5) Text summary in Command Window ----------
fprintf('\n== Summary ==\n');
for i = 1:numel(modelNames)
    fprintf('Model %-8s | logZ = %+9.3f | Posterior = %.4f\n', modelNames{i}, logZ_vec(i), post(i));
end
% fprintf('Winner: %s\n', modelNames{bestIdx});

end % main


% ===== helper: save PNG neatly =====
function savepng(figHandle, filename)
    try
        set(figHandle,'PaperPositionMode','auto');
        print(figHandle, filename, '-dpng', '-r300');
        fprintf('Saved: %s\n', filename);
    catch ME
        warning('Could not save %s: %s', filename, ME.message);
    end
end

function formatPlot()
    box on
    set(gca, 'LineWidth', 2);
    set(gca, 'FontSize', 28);
    set(gca, 'TickLabelInterpreter','latex');
    pbaspect([1.5 1 1]);
end