%% Paper benchmark formulation test
% Compares, on the paper benchmark (8 cases, target 0.4, no recalibration):
%   1. the code in main/ as is (sol_osi + gpt5_6_4_auth)
%   2. the paper formulation as written (paper_osi)
%   3. the paper formulation with the best 1, 2 and 3 changes, chosen greedily
% Score J = mean estimation error (8 cases) + mean unsafe time of cases 5-8 / 1000
% (1 % unsafe time counts like 0.001 of estimation error).
% Results: paper_benchmark_formulation_test.mat and .csv in this folder.

clear; clc
here = fileparts(mfilename('fullpath'));
cd(here)
target = 0.4;

cand = {'slow power average', struct('eta_q', 0.12); 'confirmation gate', struct('gate', true); ...
        'output smoothing', struct('smooth', true); 'lambda guard +0.02', struct('guard', 0.02); ...
        'accumulate then recover', struct('both', true); 'median-3 reference', struct('med3', true)};

rows = {'Code in main/ (as is)', runCase(struct('code', true), target); ...
        'Paper as written', runCase(struct(), target)};
chosen = struct();
left = 1:size(cand, 1);
names = {};
for k = 1:3
    best = Inf;
    for j = left
        r = runCase(mergeOpt(chosen, cand{j, 2}), target);
        fprintf('  step %d, trying %-24s J = %.4f\n', k, cand{j, 1}, r.J)
        if r.J < best, best = r.J; bj = j; br = r; end
    end
    chosen = mergeOpt(chosen, cand{bj, 2});
    left(left == bj) = [];
    names{end+1} = cand{bj, 1}; %#ok<SAGROW>
    rows(end+1, :) = {sprintf('Paper + %d: %s', k, strjoin(names, ' + ')), br}; %#ok<SAGROW>
end

get = @(f) cellfun(@(r) r.(f), rows(:, 2));
T = table(rows(:, 1), get('errW'), get('errN'), get('unsafe'), get('rough'), get('travel'), get('J'), ...
    'VariableNames', {'Method', 'ErrWideband', 'ErrNarrowband', 'UnsafeAdaptivePct', ...
    'RoughWideband', 'TravelAdaptive', 'J'});
disp(T)
save(fullfile(here, 'paper_benchmark_formulation_test.mat'), 'T')
writetable(T, fullfile(here, 'paper_benchmark_formulation_test.csv'))

function r = runCase(o, target)
evalc('R = run_paper_benchmark(target, o);');
t = R.table;
r.errW = mean(t.err([1 2 5 6]));
r.errN = mean(t.err([3 4 7 8]));
r.unsafe = mean(t.unsafe(5:8));
r.rough = mean(t.rough([1 2 5 6]));
r.travel = mean(t.travel(5:8));
r.J = mean(t.err) + r.unsafe/1000;
end

function a = mergeOpt(a, b)
for f = fieldnames(b)', a.(f{1}) = b.(f{1}); end
end