%% Paper value tuning
% Coordinate search over the seven calibrated values of the paper formulation as
% written (no structural change), on the paper benchmark, target 0.4.
% Each value is tried at x0.5, x0.7, x1.4, x2 of its current value; improvements are kept.
% Score J = mean estimation error (8 cases) + mean unsafe time of cases 5-8 / 1000.
% Results: paper_value_tuning.mat in this folder.

clear; clc
here = fileparts(mfilename('fullpath'));
cd(here)
target = 0.4;
names = {'eta_s', 'eta_A', 'l_floor', 'eta_g', 'r_g', 'g_th', 'kappa'};
cap = struct('eta_s', 0.95, 'eta_A', Inf, 'l_floor', Inf, 'eta_g', 0.95, 'r_g', 0.98, 'g_th', 6, 'kappa', Inf);
o = struct('eta_s', 0.39, 'eta_A', 4.54, 'l_floor', 1e-3, 'eta_g', 0.12, 'r_g', 0.73, 'g_th', 1.75, 'kappa', 0.625);
start = o;
[J, r0] = score(o, target);
fprintf('Table 2 values: J = %.4f\n', J)
for pass = 1:2
    for k = 1:numel(names)
        p = names{k};
        for f = [0.5 0.7 1.4 2]
            t = o;
            t.(p) = min(o.(p)*f, cap.(p));
            if strcmp(p, 'g_th'), t.g_th = max(1 + (o.g_th - 1)*f, 1.05); end
            Jt = score(t, target);
            if Jt < J - 1e-4
                J = Jt; o = t;
                fprintf('  pass %d: %-8s -> %-8.4g J = %.4f\n', pass, p, o.(p), J)
            end
        end
    end
end
[~, r1] = score(o, target);
fprintf('\n%-8s %10s %10s\n', 'value', 'Table 2', 'tuned')
for k = 1:numel(names), fprintf('%-8s %10.4g %10.4g\n', names{k}, start.(names{k}), o.(names{k})), end
fprintf('\n%-14s %8s %8s %8s %8s\n', '', 'ErrWide', 'ErrNarr', 'Unsafe%', 'RoughW')
fprintf('%-14s %8.3f %8.3f %8.1f %8.1f\n', 'Table 2', r0.errW, r0.errN, r0.unsafe, r0.rough)
fprintf('%-14s %8.3f %8.3f %8.1f %8.1f\n', 'tuned', r1.errW, r1.errN, r1.unsafe, r1.rough)
save(fullfile(here, 'paper_value_tuning.mat'), 'start', 'o', 'r0', 'r1')

function [J, r] = score(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]));
J = mean(t.err) + r.unsafe/1000;
end