%% ========================================================================
% VIBFATIGUE_REPORT_META_ACCURACY_V2
%
% PURPOSE
% -------
% This script evaluates and visualizes the accuracy of a META (adaptive)
% spectral fatigue estimator relative to:
%   (1) individual classical spectral methods, and
%   (2) a simple equal-weight averaging baseline,
% using rainflow fatigue damage as the reference.
%
% This script exists to demonstrate that discrepancies between spectral
% fatigue methods and rainflow damage are driven primarily by PSD shape
% and multimodality—not by implementation details—and to motivate the
% need for adaptive or corrected spectral estimators.
%
% The script is intended to be used *after* a Monte-Carlo PSD fatigue
% generation step and a META calibration step have been completed.
%
%
% SCIENTIFIC CONTEXT
% ------------------
% Classical spectral fatigue methods (NB, Dirlik, Zhao–Baker, Tovo–Benasciutti,
% etc.) are known to perform well only within limited PSD shape regimes
% (e.g., narrowband vs broadband, single-peak vs multimodal).
%
% This script:
%   • Quantifies method-by-method error relative to rainflow
%   • Compares META against the *best individual* method ex post
%   • Tests whether META outperforms a naive equal-weight ensemble
%   • Provides statistical evidence that adaptive weighting is required
%     when PSD shape varies
%
%
% INPUT DATA SOURCES (EXPECTED IN BASE WORKSPACE)
% -----------------------------------------------
% meta   = evalin('base','vf_meta_inputs');
%   Struct produced by the PSD Monte-Carlo fatigue generator, containing:
%     • Dspec_k4, Dspec_k6, Dspec_k8 : absolute spectral damages by method
%     • DrfOut                      : rainflow reference damages
%     • shapeOut                    : PSD shape descriptors
%     • momentsOut                  : spectral moments
%     • methodNames                 : list of spectral methods
%
% calib  = evalin('base','calib');
%   Struct produced by the META calibration script, containing:
%     • methods        : ordered list of spectral methods
%     • A              : linear weight model coefficients
%     • gate (optional): adaptive gating parameters
%     • scale          : per-k scaling factors
%     • predict_raw    : optional function handle for prediction
%
%
% DRIVER OPTIONS (optsR)
% ---------------------
% optsR.useK        : vector of fatigue exponents to include (default [4 6 8])
% optsR.makePlots   : true/false to generate figures (default true)
%
%
% FUNCTION CALL
% -------------
% results = vibfatigue_report_meta_accuracy_v2(meta, meta.DrfOut, calib, optsR);
%
%
% PRIMARY FUNCTION INPUTS
% -----------------------
% metaOrCases : struct
%   Either:
%     • metaInputs struct (preferred, from PSD MC script), or
%     • legacy cases struct array with fields:
%         .Dspec_k4, .Dspec_k6, .Dspec_k8
%
% DrfOut : [N x 4] numeric array
%   Rainflow reference damage:
%     column 1 : run index
%     column 2 : k=4 damage
%     column 3 : k=6 damage
%     column 4 : k=8 damage
%
% calib : struct
%   META calibration parameters and prediction logic
%
% opts : struct (optional)
%   Analysis and plotting options:
%     • useK            : which k values to include
%     • makePlots       : generate figures
%     • saveFigs        : save figures to disk
%     • figPrefix       : filename prefix for figures
%     • showTopN        : number of top methods to print
%     • normalizeToRF   : whether META predicts normalized damage
%     • equalAvgLabel   : label for equal-weight baseline
%
%
% OUTPUTS
% -------
% results : struct
%   Contains tables and summaries suitable for reports or publications:
%     • T4, T6, T8      : per-k error tables (all methods + META + EQAVG)
%     • Tall            : stacked results across all k
%     • S4, S6, S8      : summary statistics by method (per k)
%     • Sall            : overall summary across all k
%     • methods         : method name list
%
% Additionally, tables are optionally assigned into the base workspace:
%   vf_report_T4, vf_report_T6, vf_report_T8, vf_report_Tall, vf_report_Sall
%
%
% FIGURES GENERATED (OPTIONAL)
% ---------------------------
%   • Median absolute percent error by method
%   • CDF of |percent error| (META vs best vs EQAVG)
%   • Log-log predicted vs rainflow scatter
%   • Percent-error band plots (5/50/95%)
%
%
% INTENDED USE
% ------------
% • Course material: demonstrates why no single spectral method is sufficient
% • Handbook reference: documents a defensible META validation workflow
% • Research validation: provides statistical evidence for adaptive estimators
%
% ========================================================================


function results = vibfatigue_meta_accuracy_v2(sss, metaOrCases, DrfOut, calib, opts)
sss
%VIBFATIGUE_REPORT_META_ACCURACY_V2
% Build tables + graphs showing META accuracy vs:
%   (1) each individual method, and
%   (2) a simple equal-weight average baseline.
%
% Works with:
%   - metaInputs struct (preferred): vf_meta_inputs from your generator script
%   - legacy cases struct array with Dspec_k4/k6/k8 fields
%
% Creates:
%   - Table summaries (per k and all k stacked)
%   - Figures: boxplots, CDF/CCDF of abs % error, scatter log-log vs RF,
%              and a bar chart of median abs % error by method.
%
% Toolboxes: none (no Statistics Toolbox required). Uses basic plotting.
%
% USAGE EXAMPLES
%   meta  = evalin('base','vf_meta_inputs');
%   calib = evalin('base','calib');   % after calibration
%   results = vibfatigue_report_meta_accuracy_v2(meta, [], calib);
%
%   % or with explicit DrfOut
%   results = vibfatigue_report_meta_accuracy_v2(meta, meta.DrfOut, calib);

%% ---------------- Defaults ----------------
if nargin < 4, opts = struct; end
if nargin < 2, DrfOut = []; end

if ~isfield(opts,'useK'),        opts.useK = [4 6 8]; end
if ~isfield(opts,'makePlots'),   opts.makePlots = true; end
if ~isfield(opts,'saveFigs'),    opts.saveFigs = false; end
if ~isfield(opts,'figPrefix'),   opts.figPrefix = 'vf_meta_report'; end
if ~isfield(opts,'showTopN'),    opts.showTopN = 12; end
if ~isfield(opts,'normalizeToRF')
    if isfield(calib,'normalizeToRF'), opts.normalizeToRF = calib.normalizeToRF;
    else, opts.normalizeToRF = true;
    end
end
if ~isfield(opts,'equalAvgLabel'), opts.equalAvgLabel = 'EQAVG'; end

% If gates were trained disabled, keep them as-is in calib. Nothing to do.

%% ---------------- Ingest meta vs cases ----------------
if isstruct(metaOrCases) && numel(metaOrCases)==1 && isfield(metaOrCases,'shapeOut') && isfield(metaOrCases,'momentsOut')
    meta = metaOrCases;
    if isempty(DrfOut)
        if isfield(meta,'DrfOut'), DrfOut = meta.DrfOut; else, error('DrfOut missing.'); end
    end
    [cases, DrfOut2] = local_meta_to_cases(meta);
    if isempty(DrfOut), DrfOut = DrfOut2; end
else
    cases = metaOrCases;
    if isempty(DrfOut), error('If passing cases, you must also pass DrfOut.'); end
end

N = numel(cases);

% Determine method list from calib (best, ensures consistent ordering)
methods = calib.methods(:).';
M = numel(methods);

% Build feature matrix X (v2)
[X, runCases] = local_build_X_v2(cases);

% Align DrfOut to cases
runRf = DrfOut(:,1);
[tf, loc] = ismember(runCases, runRf);
if any(~tf), error('Some run IDs not found in DrfOut(:,1).'); end

Drf_k4 = DrfOut(loc,2);
Drf_k6 = DrfOut(loc,3);
Drf_k8 = DrfOut(loc,4);

%% ---------------- Extract spectral absolute damages (N x M) per k ----------------
Dspec4 = local_extract_Dspec_matrix(cases, 'Dspec_k4', methods);
Dspec6 = local_extract_Dspec_matrix(cases, 'Dspec_k6', methods);
Dspec8 = local_extract_Dspec_matrix(cases, 'Dspec_k8', methods);

% Equal-weight baseline (in absolute space, then optionally normalized)
wEq = ones(1,M) / M;

%% ---------------- Predict META (absolute) for each k ----------------
% If calib.normalizeToRF=true, we predict normalized yhat and multiply by RF
% Using calib.predict_raw handle if present, else local predictor.
useHandleRaw = isfield(calib,'predict_raw') && isa(calib.predict_raw,'function_handle');

yhat4 = nan(N,1); yhat6 = nan(N,1); yhat8 = nan(N,1);
yEq4  = nan(N,1); yEq6  = nan(N,1); yEq8  = nan(N,1);

for i = 1:N
    Xrow = X(i,:);
    D4 = Dspec4(i,:); D6 = Dspec6(i,:); D8 = Dspec8(i,:);
    RF4 = Drf_k4(i);  RF6 = Drf_k6(i);  RF8 = Drf_k8(i);

    % Equal-average
    yEq4(i) = sum(wEq .* D4);
    yEq6(i) = sum(wEq .* D6);
    yEq8(i) = sum(wEq .* D8);

    % META prediction
    if useHandleRaw
        yhat4(i) = calib.predict_raw(Xrow, D4, 4, RF4);
        yhat6(i) = calib.predict_raw(Xrow, D6, 6, RF6);
        yhat8(i) = calib.predict_raw(Xrow, D8, 8, RF8);
    else
        yhat4(i) = local_predict_raw(calib, Xrow, D4, 4, RF4);
        yhat6(i) = local_predict_raw(calib, Xrow, D6, 6, RF6);
        yhat8(i) = local_predict_raw(calib, Xrow, D8, 8, RF8);
    end
end

%% ---------------- Build metrics tables ----------------
% Per-k method errors + META + EQAVG
T4 = local_metrics_table(Drf_k4, Dspec4, methods, yhat4, yEq4, opts.equalAvgLabel);
T6 = local_metrics_table(Drf_k6, Dspec6, methods, yhat6, yEq6, opts.equalAvgLabel);
T8 = local_metrics_table(Drf_k8, Dspec8, methods, yhat8, yEq8, opts.equalAvgLabel);

Tall = [ ...
    local_metrics_table(Drf_k4, Dspec4, methods, yhat4, yEq4, opts.equalAvgLabel, 4); ...
    local_metrics_table(Drf_k6, Dspec6, methods, yhat6, yEq6, opts.equalAvgLabel, 6); ...
    local_metrics_table(Drf_k8, Dspec8, methods, yhat8, yEq8, opts.equalAvgLabel, 8) ];

% Summary by method
S4   = local_summarize_metrics(T4);
S6   = local_summarize_metrics(T6);
S8   = local_summarize_metrics(T8);
Sall = local_summarize_metrics(Tall);

results = struct();
results.methods = methods;
results.T4 = T4; results.T6 = T6; results.T8 = T8; results.Tall = Tall;
results.S4 = S4; results.S6 = S6; results.S8 = S8; results.Sall = Sall;

%% ---------------- Print quick comparisons ----------------
fprintf('\n=== META accuracy report (v2) ===\n');
fprintf('N=%d, methods=%d. Included k=%s\n', N, M, mat2str(opts.useK));
fprintf('Overall (all k stacked):\n');
print_top(Sall, opts.showTopN);

% Explicit improvement vs best individual + equal-average (median abs %)
medMeta = median(Sall.AbsPctMedian(strcmp(Sall.Method,'META')));
medEq   = median(Sall.AbsPctMedian(strcmp(Sall.Method,opts.equalAvgLabel)));
medBest = min(Sall.AbsPctMedian(~ismember(Sall.Method, {'META',opts.equalAvgLabel})));

fprintf('\nMedian abs %% error:\n');
fprintf('  META   : %.3f%%\n', medMeta);
fprintf('  EQAVG  : %.3f%%\n', medEq);
fprintf('  BestInd: %.3f%%\n', medBest);

fprintf('Improvement (lower is better):\n');
fprintf('  META vs EQAVG  : %.2fx\n', medEq/max(realmin,medMeta));
fprintf('  META vs BestInd: %.2fx\n', medBest/max(realmin,medMeta));

%% ---------------- Plots ----------------
if opts.makePlots
    % 1) Bar chart of median abs % error (all k)
    figure('Name','Median Abs % Error (All k)','Color','w');
    SallPlot = Sall;
    % order by median abs % error
    [~,ord] = sort(SallPlot.AbsPctMedian,'ascend');
    SallPlot = SallPlot(ord,:);
    bar(SallPlot.AbsPctMedian);
    set(gca,'XTick',1:height(SallPlot),'XTickLabel',SallPlot.Method,'XTickLabelRotation',45);
    ylabel('Median |% error|');
    title('Median absolute percent error vs RF (All k stacked)');
    grid on;

    % 2) CDF of abs % error for META vs best individual vs EQAVG (all k)
    [bestName, ~] = best_individual(Sall);
    figure('Name','CDF Abs % Error','Color','w');
    hold on;
    cdf_plot(abs(Tall.PctErr(strcmp(Tall.Method,'META'))));
    cdf_plot(abs(Tall.PctErr(strcmp(Tall.Method,opts.equalAvgLabel))));
    cdf_plot(abs(Tall.PctErr(strcmp(Tall.Method,bestName))));
    hold off;
    legend({'META',opts.equalAvgLabel,bestName},'Location','southeast');
    xlabel('|% error|'); ylabel('CDF');
    title('CDF of absolute percent error (All k stacked)');
    grid on;

    % 3) Log-log scatter: predicted vs RF (META + best + EQAVG)
    figure('Name','Predicted vs RF (All k)','Color','w');
    hold on;
    scatter_loglog(Tall.RF(strcmp(Tall.Method,'META')), Tall.Pred(strcmp(Tall.Method,'META')));
    scatter_loglog(Tall.RF(strcmp(Tall.Method,opts.equalAvgLabel)), Tall.Pred(strcmp(Tall.Method,opts.equalAvgLabel)));
    scatter_loglog(Tall.RF(strcmp(Tall.Method,bestName)), Tall.Pred(strcmp(Tall.Method,bestName)));
    plot_identity_line(Tall.RF);
    hold off;
    legend({'META',opts.equalAvgLabel,bestName,'y=x'},'Location','best');
    xlabel('RF damage'); ylabel('Predicted damage');
    title('Predicted vs RF (log-log, all k stacked)');
    grid on;

    % 4) Boxplot-like summary without toolbox: use quantile bands per method (all k)
    figure('Name','Error Bands (All k)','Color','w');
    local_error_band_plot(Tall, methods, opts.equalAvgLabel);
    title('Percent error bands (5/50/95%) by method (All k stacked)');
    ylabel('% error'); grid on;

    if opts.saveFigs
        saveas(gcf, [opts.figPrefix '_errbands.png']);
    end
end

%% ---------------- Return also as workspace tables if desired ----------------
try
    assignin('base','vf_report_T4',T4);
    assignin('base','vf_report_T6',T6);
    assignin('base','vf_report_T8',T8);
    assignin('base','vf_report_Tall',Tall);
    assignin('base','vf_report_Sall',Sall);
    fprintf('\nAssigned to base: vf_report_T4/T6/T8/Tall and vf_report_Sall\n');
catch
end

sss

writetable(Sall, sss);

end

%% ========================== helpers ==========================

function [cases, DrfOut] = local_meta_to_cases(meta)
shape = meta.shapeOut;
mom   = meta.momentsOut;

run = shape(:,1);
nu0 = shape(:,2);
nu1 = shape(:,3);
alpha = shape(:,4);
epsBW = shape(:,5);
SSI = shape(:,6);
SCD = shape(:,7);
Hn = shape(:,9);
domFrac = shape(:,11);
Neff = shape(:,12);

m0 = mom(:,2); m1 = mom(:,3); m2 = mom(:,4); m4 = mom(:,5);

methods = meta.methodNames(:);
N = numel(run);
M = numel(methods);

D4 = meta.Dspec_k4; D6 = meta.Dspec_k6; D8 = meta.Dspec_k8;

cases(N,1) = struct();
for i = 1:N
    cases(i).run = run(i);
    cases(i).nu0 = nu0(i);
    cases(i).nu1 = nu1(i);
    cases(i).alpha = alpha(i);
    cases(i).epsBW = epsBW(i);
    cases(i).SSI = SSI(i);
    cases(i).SCD = SCD(i);
    cases(i).Hn = Hn(i);
    cases(i).domFrac = domFrac(i);
    cases(i).Neff = Neff(i);

    cases(i).m0 = m0(i);
    cases(i).m1 = m1(i);
    cases(i).m2 = m2(i);
    cases(i).m4 = m4(i);

    s4 = struct(); s6 = struct(); s8 = struct();
    for j = 1:M
        nm = methods{j};
        s4.(nm) = D4(i,j);
        s6.(nm) = D6(i,j);
        s8.(nm) = D8(i,j);
    end
    cases(i).Dspec_k4 = s4;
    cases(i).Dspec_k6 = s6;
    cases(i).Dspec_k8 = s8;
end

DrfOut = meta.DrfOut;
end

function [X, runCases] = local_build_X_v2(cases)
N = numel(cases);
runCases = (1:N).';
if isfield(cases(1),'run')
    runCases = arrayfun(@(c)c.run, cases(:));
end

% Required / optional fields
epsBW   = get_or_nan(cases,'epsBW',N);
Hn      = get_or_nan(cases,'Hn',N);
Neff    = get_or_nan(cases,'Neff',N);
alpha   = get_or_nan(cases,'alpha',N);
domFrac = get_or_nan(cases,'domFrac',N);
SSI     = get_or_nan(cases,'SSI',N);
SCD     = get_or_nan(cases,'SCD',N);

hasMom = isfield(cases(1),'m0') && isfield(cases(1),'m1') && isfield(cases(1),'m2') && isfield(cases(1),'m4');

% Fill alpha if needed from moments
if any(~isfinite(alpha)) && hasMom
    m0 = arrayfun(@(c)c.m0, cases(:));
    m2 = arrayfun(@(c)c.m2, cases(:));
    m4 = arrayfun(@(c)c.m4, cases(:));
    nu0 = sqrt(m2 ./ max(m0,realmin));
    nu1 = sqrt(m4 ./ max(m2,realmin));
    aHat = nu0 ./ max(nu1,realmin);
    idx = ~isfinite(alpha);
    alpha(idx) = aHat(idx);
end

% epsBW from alpha where missing
if any(~isfinite(epsBW))
    idx = ~isfinite(epsBW);
    epsBW(idx) = sqrt(max(0, 1 - alpha(idx).^2));
end

% SSI/SCD from moments if needed
if (any(~isfinite(SSI)) || any(~isfinite(SCD))) && hasMom
    m0 = arrayfun(@(c)c.m0, cases(:));
    m2 = arrayfun(@(c)c.m2, cases(:));
    m4 = arrayfun(@(c)c.m4, cases(:));
    SSIc = (m0 .* m4) ./ max(m2.^2, realmin);
    SSI(~isfinite(SSI)) = SSIc(~isfinite(SSI));
    SCDc = sqrt(max(0, SSIc));
    SCD(~isfinite(SCD)) = SCDc(~isfinite(SCD));
end

alpha1 = nan(N,1); alpha2 = nan(N,1); delta = nan(N,1);
if hasMom
    m0 = arrayfun(@(c)c.m0, cases(:));
    m1 = arrayfun(@(c)c.m1, cases(:));
    m2 = arrayfun(@(c)c.m2, cases(:));
    m4 = arrayfun(@(c)c.m4, cases(:));
    alpha1 = m1 ./ sqrt(max(m0.*m2, realmin));
    alpha2 = m2 ./ sqrt(max(m0.*m4, realmin));
    delta  = sqrt(max(0, 1 - alpha1));
end

Neff01 = clamp01((Neff - 1)/4);

X = [ones(N,1), epsBW(:), Hn(:), Neff01(:), alpha(:), domFrac(:), SSI(:), SCD(:), alpha1(:), alpha2(:), delta(:)];
end

function v = get_or_nan(cases, field, N)
if isfield(cases(1), field)
    v = arrayfun(@(c)c.(field), cases(:));
else
    v = nan(N,1);
end
end

function D = local_extract_Dspec_matrix(cases, specField, methods)
N = numel(cases);
M = numel(methods);
D = nan(N,M);
for m = 1:M
    nm = methods{m};
    if ~isfield(cases(1).(specField), nm)
        error('Missing field cases(1).%s.%s', specField, nm);
    end
    D(:,m) = arrayfun(@(c)c.(specField).(nm), cases(:));
end
end

function yhat_raw = local_predict_raw(calib, Xrow, DspecRow, kVal, RF)
Dnorm = DspecRow ./ max(RF, calib.epsNorm);
yhat_norm = local_predict_norm(calib, Xrow, Dnorm, kVal);
yhat_raw = yhat_norm .* RF;
end

function yhat = local_predict_norm(calib, Xrow, DnormRow, kVal)
X = Xrow(:).';
D = DnormRow(:).';

S = (calib.A * X')';
w = softmax_rows(S);

% v2: alpha col 5, domFrac col 6
alpha = X(5);
domFrac = X(6);

if isfield(calib,'gate') && isfield(calib.gate,'enabled') && calib.gate.enabled
    gA = calib.gate.gMaxA * sigmoid(calib.gate.sA * (alpha   - calib.gate.thA));
    gD = calib.gate.gMaxD * sigmoid(calib.gate.sD * (domFrac - calib.gate.thD));
    gi = calib.gate.gateIdx;
    w(gi) = w(gi) + gA + gD;
    w = w / sum(w);
end

mix = sum(w(:)'.*D);

Sk = 1;
if isfield(calib,'scale')
    if kVal==4, Sk = calib.scale.S4; end
    if kVal==6, Sk = calib.scale.S6; end
    if kVal==8, Sk = calib.scale.S8; end
end
yhat = Sk * mix;
end

function T = local_metrics_table(RF, DspecMat, methods, yMETA, yEQ, eqLabel, kValOpt)
if nargin < 7, kValOpt = []; end
N = numel(RF);
M = numel(methods);

% Stack rows: each method, plus META and EQAVG
nRows = N*(M+2);

Method = strings(nRows,1);
RFv    = nan(nRows,1);
Pred   = nan(nRows,1);
PctErr = nan(nRows,1);

row = 0;

% Individual methods
for m = 1:M
    y = DspecMat(:,m);
    for i = 1:N
        row = row + 1;
        Method(row) = string(methods{m});
        RFv(row) = RF(i);
        Pred(row) = y(i);
        PctErr(row) = 100*(y(i)/max(realmin,RF(i)) - 1);
    end
end

% EQAVG
for i = 1:N
    row = row + 1;
    Method(row) = string(eqLabel);
    RFv(row) = RF(i);
    Pred(row) = yEQ(i);
    PctErr(row) = 100*(yEQ(i)/max(realmin,RF(i)) - 1);
end

% META
for i = 1:N
    row = row + 1;
    Method(row) = "META";
    RFv(row) = RF(i);
    Pred(row) = yMETA(i);
    PctErr(row) = 100*(yMETA(i)/max(realmin,RF(i)) - 1);
end

T = table(Method, RFv, Pred, PctErr, 'VariableNames', {'Method','RF','Pred','PctErr'});
if ~isempty(kValOpt)
    T.k = repmat(kValOpt, height(T), 1);
end
end

function S = local_summarize_metrics(T)
methods = unique(T.Method,'stable');
m = numel(methods);

Method = strings(m,1);
N = zeros(m,1);
AbsPctMedian = nan(m,1);
AbsPctMean   = nan(m,1);
Pct05 = nan(m,1);
Pct95 = nan(m,1);
RMSElog = nan(m,1);

for i = 1:m
    Method(i) = methods(i);
    idx = T.Method==methods(i);
    e = T.PctErr(idx);
    e = e(isfinite(e));
    N(i) = numel(e);
    if N(i)==0, continue; end

    AbsPctMedian(i) = median(abs(e));
    AbsPctMean(i)   = mean(abs(e));
    Pct05(i) = prctile_basic(e,5);
    Pct95(i) = prctile_basic(e,95);

    % rmse in log space (Pred vs RF)
    RF = T.RF(idx);  P = T.Pred(idx);
    good = isfinite(RF) & isfinite(P) & RF>0 & P>0;
    r = log(P(good)) - log(RF(good));
    RMSElog(i) = sqrt(mean(r.^2));
end

S = table(Method, N, AbsPctMedian, AbsPctMean, Pct05, Pct95, RMSElog);
end

function print_top(S, topN)
[~,ord] = sort(S.AbsPctMedian,'ascend');
topN = min(topN, numel(ord));
for i = 1:topN
    j = ord(i);
    fprintf('  %-6s  med|%%|=%7.3f  mean|%%|=%7.3f  RMSElog=%7.4f  (N=%d)\n', ...
        S.Method(j), S.AbsPctMedian(j), S.AbsPctMean(j), S.RMSElog(j), S.N(j));
end
end

function [bestName, bestVal] = best_individual(Sall)
mask = ~ismember(Sall.Method, ["META","EQAVG"]);
vals = Sall.AbsPctMedian(mask);
names = Sall.Method(mask);
[bestVal,ix] = min(vals);
bestName = char(names(ix));
end

function cdf_plot(x)
x = x(isfinite(x) & x>=0);
x = sort(x);
n = numel(x);
if n==0, return; end
y = (1:n)'/n;
plot(x,y,'LineWidth',1.5);
end

function scatter_loglog(x,y)
x = x(:); y = y(:);
good = isfinite(x) & isfinite(y) & x>0 & y>0;
loglog(x(good), y(good), '.', 'MarkerSize', 8);
end

function plot_identity_line(x)
x = x(isfinite(x) & x>0);
if isempty(x), return; end
xmin = min(x); xmax = max(x);
loglog([xmin xmax],[xmin xmax],'-','LineWidth',1.0);
end

function local_error_band_plot(Tall, methods, eqLabel)
% Build list = individual methods + EQAVG + META
names = [methods(:); {eqLabel}; {'META'}];
n = numel(names);

% y positions
hold on;
for i = 1:n
    nm = string(names{i});
    e = Tall.PctErr(Tall.Method==nm);
    e = e(isfinite(e));
    if isempty(e), continue; end

    q5  = prctile_basic(e,5);
    q50 = prctile_basic(e,50);
    q95 = prctile_basic(e,95);

    y0 = n - i + 1;

    % whisker line (q5 to q95)
    plot([q5 q95],[y0 y0],'-','LineWidth',2);
    % median marker
    plot(q50, y0, 'o', 'MarkerSize', 6, 'LineWidth', 1.5);
end
set(gca,'YTick',1:n,'YTickLabel',flip(string(names)),'YDir','normal');
xlabel('% error');
end

function p = prctile_basic(x, q)
x = x(isfinite(x));
x = sort(x(:));
n = numel(x);
if n==0, p = NaN; return; end
if q<=0, p = x(1); return; end
if q>=100, p = x(end); return; end
r = 1 + (n-1)*(q/100);
lo = floor(r); hi = ceil(r);
if lo==hi
    p = x(lo);
else
    p = x(lo) + (r-lo)*(x(hi)-x(lo));
end
end

function W = softmax_rows(S)
S = S - max(S,[],2);
E = exp(S);
W = E ./ sum(E,2);
end

function y = sigmoid(x)
y = 1 ./ (1 + exp(-x));
end

function y = clamp01(x)
y = min(1, max(0, x));
end
