function R = run_code_bridge_benchmark(target, diag, seed_offset, idle)
%RUN_PAPER_BENCHMARK Paper 2 benchmark: the eight calibration cases of Table 3.
%   R = run_paper_benchmark        % target margin 0.4 (Paper 1 guideline)
%   R = run_paper_benchmark(0.55)
%
% Closed loop of the discrete outer PID and the two identified plants. Damage
% is imposed by mixing the plant outputs, y = (1-delta) y0 + delta yd, so the
% true NCF margin of lambda*C*[(1-delta)H0 + delta Hd] is known at every update.
% Estimator and adaptive law: paper_osi, the equations of Section 3.1 with Table 2.
% Metrics (t >= 5 s):
%   err   mean |b_est - b_true|                     estimation accuracy
%   rough 1000*median|b(n) - 2b(n-1) + b(n-2)|      estimate smoothness
%   travel sum|lambda(n) - lambda(n-1)|             controller smoothness
%   unsafe percent of time with b_true < target - 0.02 (active cases)

if nargin < 4, idle = false; end
if nargin < 3, seed_offset = 0; end
if nargin < 1, target = 0.4; end
if nargin < 2, diag = struct(); end
here = fileparts(mfilename('fullpath'));
main = fileparts(here);

dt = 1/2048;
E0 = load(fullfile(main, 'saved_results', 'exp_tf', 'tf_sp1_7000_20260627_210411.mat'));
Ed = load(fullfile(main, 'saved_results', 'exp_tf', 'tf_sp1_5000_20260627_233333.mat'));
pv = Ed.pid_values;
C = pid(pv.P, pv.I, pv.D, 1/pv.derivative_filter_N, dt, ...
    'IFormula', 'ForwardEuler', 'DFormula', 'ForwardEuler');
H0 = ss(c2d(E0.H_dc_dm, dt, 'zoh'));
Hd = ss(c2d(Ed.H_dc_dm, dt, 'zoh'));
fg = (2:0.25:80)';
fg = fg(fg < 58 | fg > 62);
L0f = squeeze(freqresp(C*H0, 2*pi*fg));
Ldf = squeeze(freqresp(C*Hd, 2*pi*fg));

dur = 120;
hist = {'Sudden', 'Gradual'};
bands = {[2 80], [2 10]};
bandName = {'Wideband', 'Narrowband'};
k = 0;
for adapt = [false true]
    for ib = 1:2
        for ih = 1:2
            k = k + 1;
            rng(seed_offset + 10*ib + ih, 'twister')   % frozen and active cases share signals
            cs = struct('id', k, 'history', hist{ih}, 'band', bandName{ib}, ...
                'adaptive', adapt);
            R.cases(k) = simulate(cs, bands{ib}, dur, dt, C, H0, Hd, L0f, Ldf, target, diag, idle);
        end
    end
end
R.target = target;
R.table = struct2table(rmfield(R.cases, {'t', 'dc', 'dm', 'tc', 'delta', 'b_est', 'b_true', 'lambda'}));
disp(R.table)
R.seed_offset = seed_offset; R.idle = idle;
end

function c = simulate(c, band, dur, dt, C, H0, Hd, L0f, Ldf, target, diag, idle)
N = round(dur/dt);
t = (0:N-1)'*dt;
if c.history(1) == 'S'
    delta = double(t >= 40 & t < 80);
else
    delta = interp1([0 20 50 70 100 dur], [0 0 1 1 0 0], t);
end
u = bandNoise(N, 1/dt, band)*0.15;
v = 0.05*randn(N, 1);
if idle, u(t>=55 & t<70)=0; end

prm = struct('fs', 1/dt, 'dt', 0.125, 'Tw', 0.5, 'band', [2 80], 'eta_s', 0.39, ...
    'eta_A', 4.54, 'sigma_u', 0.0575/4.54, 'l_floor', 1e-3, 'eta_g', 0.12, 'r_g', 0.73, ...
    'g_th', 1.75, 'b_tar', target, 'kappa', 0.625, 'lambda_min', 0.2, 't_ref', 1.5);
if isfield(diag, 'code')
    addpath(fullfile(fileparts(fileparts(mfilename('fullpath'))), 'control_and_HS', 'OSI'))
    scfg = sol_osi('config', dt, C*H0, struct()); sst = sol_osi('state', scfg);
    aprm = gpt5_6_4_auth('params', struct('target', target, 'dt_update', 0.125)); ast = gpt5_6_4_auth('state', aprm);
elseif isfield(diag, 'sol')
    scfg = sol_osi_ablate('config', dt, C*H0, diag.sol); sst = sol_osi_ablate('state', scfg);
else
    for fn = fieldnames(diag)', prm.(fn{1}) = diag.(fn{1}); end
end
ost = paper_osi('init', @(f) squeeze(freqresp(C*H0, 2*pi*f)), prm);

[Ac, Bc, Cc, Dc] = ssdata(ss(C));
[A0, B0, C0] = ssdata(H0);
[Ad, Bd, Cd] = ssdata(Hd);
xc = zeros(size(Ac, 1), 1); x0 = zeros(size(A0, 1), 1); xd = zeros(size(Ad, 1), 1);
nb = ost.nb;
nU = floor(N/nb);
[tc, bE, bT, lam] = deal(NaN(nU, 1));
dc = zeros(N, 1); dm = zeros(N, 1);
a = 1;
for i = 1:N
    dm(i) = (1 - delta(i))*(C0*x0) + delta(i)*(Cd*xd) + v(i);
    e = u(i) - dm(i);
    dc(i) = a*(Cc*xc + Dc*e);
    xc = Ac*xc + Bc*e;
    x0 = A0*x0 + B0*dc(i);
    xd = Ad*xd + Bd*dc(i);
    if mod(i, nb) == 0
        n = i/nb;
        if isfield(diag, 'code')                                    % code in main/ as is
            sst.a_cmd = a; [b, sst, oi] = sol_osi(u(i-nb+1:i), dm(i-nb+1:i), scfg, sst);
            if c.adaptive, [ost.lambda, ast] = gpt5_6_4_auth(b, oi.full_margin, oi.residual_weight, ast, aprm); end
        elseif isfield(diag, 'sol')                                 % code estimator (diagnosis)
            sst.a_cmd = 1; [b, sst] = sol_osi_ablate(u(i-nb+1:i), dm(i-nb+1:i), scfg, sst);
        else
            [b, ost] = paper_osi(u(i-nb+1:i), dm(i-nb+1:i), ost);
        end
        if ~c.adaptive, ost.lambda = 1; end
        La = a*((1 - delta(i))*L0f + delta(i)*Ldf);
        tc(n) = t(i); bE(n) = b; lam(n) = a;
        a = ost.lambda;
        bT(n) = min(abs(1 + La)./sqrt(2*(1 + abs(La).^2)));
    end
end
w = tc >= 5;
w = w & ~isnan(bE);
c.err = mean(abs(bE(w) - bT(w)));
c.rough = 1000*median(abs(diff(bE(w), 2)));
c.travel = sum(abs(diff(lam(w))));
c.unsafe = 100*mean(bT(w) < target - 0.02);
s = 1:8:N;   % keep 256 Hz for plotting
c.t = t(s); c.dc = dc(s); c.dm = dm(s);
c.tc = tc; c.delta = interp1(t, delta, tc); c.b_est = bE; c.b_true = bT; c.lambda = lam;
end

function x = bandNoise(N, fs, band)
X = fft(randn(N, 1));
f = (0:N-1)'*fs/N;
f = min(f, fs - f);
x = real(ifft(X.*(f >= band(1) & f <= band(2))));
x = x/rms(x);
end