function [b, st] = paper_osi(u, y, st)
%PAPER_OSI NCF OSI and adaptive factor of Paper 2 Section 3.1 (defaults: the first draft).
%   Current text: c_ref, c_gate, c_paper, up_gain 0.32, d_mode 'c' (see run_paper_benchmark).
%   st = paper_osi('init', L0fun, prm)   L0fun(f_hz) returns the initial loop C*H0
%   [b, st] = paper_osi(u_blk, y_blk, st) consumes one update block (Delta t)
%   and returns b_OSI; st.lambda holds lambda_{n+1}.
% Only deviation from the text: the second-smallest bin of beta is reported.

if ischar(u)
    p = st;                                    % parameters (Table 2)
    N = p.Tw*p.fs;
    f = (0:N/2)'*p.fs/N;
    k = f >= p.band(1) & f <= p.band(2) & ~(f >= 58 & f <= 62);
    st = p;
    st.N = N; st.nb = p.dt*p.fs; st.bins = find(k); st.f = f(k);
    st.w = 0.5*(1 - cos(2*pi*(0:N-1)'/(N - 1)));                    % Hann
    st.L0 = y(st.f);
    st.G = logspace(log10(0.5), log10(6), 25);
    st.ub = zeros(N, 1); st.yb = zeros(N, 1); st.filled = 0; st.n = 0;
    st.puu = []; st.pyy = []; st.pyu = [];
    st.pyy0 = 0; st.nref = 0; st.lgs = zeros(size(st.f));
    st.lambda = 1;
    if ~isfield(st, 'recovery'), st.recovery = 'median'; end   % diagnosis switches
    if ~isfield(st, 'eta_q'), st.eta_q = st.eta_s; end
    if ~isfield(st, 'gfit'), st.gfit = 'bin'; end
    if ~isfield(st, 'gate'), st.gate = false; end
    if ~isfield(st, 'both'), st.both = false; end
    if ~isfield(st, 'med3'), st.med3 = false; end
    if ~isfield(st, 'smooth'), st.smooth = false; end
    if ~isfield(st, 'guard'), st.guard = 0; end
    if ~isfield(st, 'drive'), st.drive = false; end       % replace power branch by drive ratio
    if ~isfield(st, 'drive_avg'), st.drive_avg = false; end
    if ~isfield(st, 'quiet_fit'), st.quiet_fit = false; end
    if ~isfield(st, 'drive_floor'), st.drive_floor = false; end
    if ~isfield(st, 'deadband'), st.deadband = 0; end
    if ~isfield(st, 'up_gain'), st.up_gain = 1; end        % relative upward lambda gain
    if ~isfield(st, 'c_paper'), st.c_paper = false; end    % v109: c_n = (rho-1)/(rho_th-1) in [0.1,1]
    if ~isfield(st, 'd_mode'), st.d_mode = 'ramp'; end     % 'ramp' (d_n), 'one' (no d_n), 'c' (d_n = c_n)
    % Atomic ports from sol_osi / gpt5_6_4_auth, all off by default.
    flags = {'c_norm','c_gate','c_level','c_prior','c_den','c_idmed', ...
        'c_powmed','c_ref','c_activity','c_kick','c_cap','c_guard', ...
        'c_deadband','c_up','c_rate','c_step'};
    for j = 1:numel(flags)
        if ~isfield(st,flags{j}), st.(flags{j}) = false; end
    end
    st.rho_drive = 1; st.ratio_prev = 0;
    st.bsm = [];
    st.snap = zeros(numel(st.f), 3);
    st.pyyq = [];
    b = st;
    return
end

nb = numel(u);
st.ub = [st.ub(nb+1:end); u(:)];
st.yb = [st.yb(nb+1:end); y(:)];
st.filled = min(st.filled + nb, st.N);
b = NaN;
if st.filled < st.N, return, end
st.n = st.n + 1;

U = fft(st.w.*st.ub); U = U(st.bins);
Y = fft(st.w.*st.yb); Y = Y(st.bins);
puu = abs(U).^2; pyy = abs(Y).^2; pyu = Y.*conj(U);
if isempty(st.puu)
    st.puu = puu; st.pyy = pyy; st.pyu = pyu;
else
    e = st.eta_s;
    st.puu = (1 - e)*st.puu + e*puu;
    st.pyy = (1 - e)*st.pyy + e*pyy;
    st.pyu = (1 - e)*st.pyu + e*pyu;
end

% direct identification
T = st.pyu./st.puu;
g2 = abs(st.pyu).^2./(st.puu.*st.pyy);
L = T./(1 - T);
Ls = complex(movmedian(real(L), 3), movmedian(imag(L), 3));
bid = chordal(Ls);
A = rms(st.ub) > st.eta_A*st.sigma_u;
lraw = st.puu/max(st.puu);
wid = A*g2.*lraw./(lraw + st.l_floor);
if st.c_prior || st.c_den
    Tu = T;
    if st.c_prior
        Lp = st.lambda*st.L0;
        Tu = wid.*T + (1-wid).*Lp./(1+Lp);
    end
    den = 1-Tu;
    if st.c_den
        small = abs(den)<0.02;
        den(small) = 0.02*exp(1i*angle(den(small)+eps));
    end
    L = Tu./den;
    Ls = complex(movmedian(real(L),3),movmedian(imag(L),3));
    bid = chordal(Ls);
end
if st.c_idmed, bid = movmedian(chordal(L),3); end
drive_ratio = 1;
if st.c_gate || st.c_kick
    w2 = wid.^2; Lm = st.lambda*st.L0;
    rho_code = max(real(sum(w2.*conj(Lm).*L))/max(sum(w2.*abs(Lm).^2),eps),eps);
    st.cn = min(max(2*(rho_code-1),0.1),1);
end
if st.drive || st.drive_avg || st.drive_floor || strcmp(st.recovery, 'drive')
    Lm = st.lambda*st.L0; w2 = wid.^2;
    rho = real(sum(w2.*conj(Lm).*L))/max(sum(w2.*abs(Lm).^2), eps);
    rho = min(max(rho, st.G(1)), st.G(end));
end

% output-power inference (reference: first t_ref seconds at lambda = 1)
if st.drive || st.drive_avg
    if st.n*st.dt > st.t_ref
        if st.drive
            st.lgs(:) = log(rho);
        else
            st.lgs = (1 - st.eta_g)*st.lgs + st.eta_g*log(rho);
        end
    end
else
if isempty(st.pyyq), st.pyyq = pyy; else, st.pyyq = (1 - st.eta_q)*st.pyyq + st.eta_q*pyy; end
if st.n*st.dt <= st.t_ref
    st.nref = st.nref + 1;
    if st.c_ref                                                     % code snapshots at t=0.5,1.0,1.5 s
        step = st.n + st.N/st.nb - 1;
        if mod(step,4)==0 && step<=12, st.snap(:,step/4) = pyy; end
        st.pyy0 = median(st.snap,2);
    elseif st.med3
        if mod(st.nref, 4) == 0 && st.nref <= 12, st.snap(:, st.nref/4) = pyy; end
        st.pyy0 = median(st.snap, 2);
    else
        st.pyy0 = st.pyy0 + (st.pyyq - st.pyy0)/st.nref;
    end
else
    qobs = st.pyyq./st.pyy0;
    reference_ratio = median(qobs);
    if st.c_norm, qobs = qobs/max(reference_ratio^0.5,eps); end
    if st.c_kick
        raw_kick = max(rho_code/st.rho_drive-1,0);
        kick = min(min(raw_kick,st.ratio_prev),0.95);
        st.ratio_prev = raw_kick;
        drive_ratio = 1+0.8*A*kick;
        st.rho_drive = st.rho_drive+(1-exp(-st.dt/2))*(rho_code-st.rho_drive);
    end
    switch st.recovery
        case 'median', rec = median(qobs) < 1;
        case 'none', rec = false;
        case 'drive', rec = rho <= 1;
    end
    if rec && ~st.both
        st.lgs = st.r_g*st.lgs;                                     % recovery
    else
        qmod = abs((1 + st.L0)./(1 + st.lambda*st.L0.*st.G)).^2;
        if st.c_norm, qmod = qmod./max(median(qmod,1).^0.5,eps); end
        if strcmp(st.gfit, 'joint')                                  % one g for all bins
            [~, i] = min(sum(abs(log(qobs) - log(qmod)), 1));
            i = repmat(i, size(qobs));
        else
            [~, i] = min(abs(log(qobs) - log(qmod)), [], 2);
        end
        lg = log(st.G(i)');
        if st.c_gate || st.c_level
            c = 1;
            if st.c_gate, c = min(max(2*(rho_code-1),0),1)+0.1*A*(1-min(max(2*(rho_code-1),0),1)); end
            if st.c_paper, c = st.cn; end
            if st.c_level && st.c_gate
                % Same extra level-dependent fit allowance as the code.
                coherent = 0;
                if st.c_gate, coherent = min(max(2*(rho_code-1),0),1); end
                c = c+0.28*max(1-4/max(reference_ratio,eps),0)*A*(1-coherent);
            end
            lg = log(1+c*(exp(lg)-1));
        end
        if st.quiet_fit, lg(wid >= 0.5) = 0; end
        if st.gate                                                    % drive-ratio gate
            w2 = wid.^2; Lm = st.lambda*st.L0;
            rho = real(sum(w2.*conj(Lm).*L))/max(sum(w2.*abs(Lm).^2), eps);
            c = min(max(2*(rho - 1), 0), 1);
            lg = (c + 0.1*(1 - c))*lg;
        end
        eg = st.eta_g;
        if st.c_activity, eg = A*eg; end
        st.lgs = (1 - eg)*st.lgs + eg*lg;
        if rec && st.both, st.lgs = st.r_g*st.lgs; end             % code order: accumulate, then recover
    end
end
end
if st.drive_floor && st.n*st.dt > st.t_ref
    st.lgs = max(st.lgs, log(rho));
end
gs = exp(st.lgs);
gs = gs*drive_ratio;
bpow = chordal(st.lambda*gs.*st.L0);
if st.c_powmed, bpow = movmedian(bpow,3); end
d = min(max((gs - 1)/(st.g_th - 1), 0), 1);
switch st.d_mode
    case 'one', d = ones(size(gs));
    case 'c', d = st.cn*ones(size(gs));
end
wpow = (1 - wid).*d;

% prior and per-frequency fusion
b0 = chordal(st.lambda*st.L0);
if st.c_cap
    bid = min(bid,1/sqrt(2)); bpow = min(bpow,1/sqrt(2)); b0 = min(b0,1/sqrt(2));
end
w0 = (1 - wid).*(1 - d);
beta = wid.*bid + wpow.*bpow + w0.*b0;
bs = sort(beta);
[~, order] = sort(beta);
st.residual_weight = wpow(order(2));
st.activity = A;
st.g_memory = exp(st.lgs);
b = bs(2);
if st.smooth                                                        % output smoothing, tau = 0.25 s
    if isempty(st.bsm), st.bsm = b; end
    st.bsm = st.bsm + (1 - exp(-st.dt/0.25))*(b - st.bsm); b = st.bsm;
end

% adaptive factor: each code detail is independently switchable.
if st.c_guard || st.c_deadband || st.c_up || st.c_rate || st.c_step
    a0 = st.lambda;
    b_ref = st.b_tar+st.guard;
    if st.c_guard, b_ref = st.b_tar+0.02*(1+st.residual_weight); end
    e = b-b_ref;
    db = st.deadband;
    if st.c_deadband, db = 0.008*(1+st.residual_weight); end
    e = sign(e)*max(abs(e)-db,0);
    up = st.up_gain;
    if st.c_up, up = 2*b_ref^2; end
    gain = st.kappa;
    if st.c_rate, gain = gain*sqrt(a0); end
    a_free = min(max(a0*exp(gain*(min(e,0)+up*max(e,0))),st.lambda_min),1);
    st.lambda = a_free;
    if st.c_step, st.lambda = a0+min(max(a_free-a0,-0.05),0.05); end
elseif st.deadband == 0 && st.up_gain == 1
    st.lambda = min(max(st.lambda*exp(st.kappa*(b - st.b_tar - st.guard)), st.lambda_min), 1);
else
    e = b - st.b_tar - st.guard;
    e = sign(e)*max(abs(e) - st.deadband, 0);
    e = min(e, 0) + st.up_gain*max(e, 0);
    st.lambda = min(max(st.lambda*exp(st.kappa*e), st.lambda_min), 1);
end
end

function c = chordal(L)
c = abs(1 + L)./sqrt(2*(1 + abs(L).^2));
end