function [Wbar_imp, Ybar_imp, Nbar_imp, Wof_after_all, Zf_after_all, Yf_after_all] = ge_forward_eval_aftershock( ...
    xvec, params, alphaF, alphaOf_all, alphaJ_all, Q_after_all, W_before_all, zbar0, Ybar0, Wof_fun)

% Implied aggregates given (log gW, log gY); parallel across firms.

zeta  = params.zeta; chi = params.chi; theta = params.theta;
phi   = params.phi;  sigma = params.sigma; Phi = params.Phi; 

gW = exp(xvec(1));
gY = exp(xvec(2));
gN = gW^(phi*(1-sigma)/(1+phi*sigma));

% zbar multiplier (“blue” factor)
S  = gW^(1 + phi - zeta) * gN^(-phi*sigma) / gY;

F = numel(alphaF);
Wof_after_all = cell(F,1);
Zf_after_all  = nan(F,1);

% ---------- parallel per-firm solves ----------
parfor f = 1:F
    aOf = alphaOf_all{f};
    if isempty(aOf), continue; end

    params_f           = params;
    params_f.zbar      = zbar0 * S;   % scale LHS
    params_f.rhs_scale = gY;          % scale RHS (patch in solve_task_wages)

    w0_after = W_before_all{f};       % warm start

    [W_after, Zf_after, Z_component_after, details_after] = solve_task_wages(Q_after_all{f}, alphaJ_all{f}, aOf, alphaF(f),params_f, w0_after);
    Zf_after_all(f) = Zf_after;

    % NEW: collect primal Y_f
    if isfield(details_after,'Yf') && ~isempty(details_after.Yf)
        Yf_after_all(f) = details_after.Yf;
    else
        % (very unlikely if patch #2 is in place) recompute on the fly
        [Yf_after_all(f), ~] = Yf_from_wq(W_after, Q_after_all{f}, alphaJ_all{f}, alphaOf_all{f}, alphaF(f),params);
    end

    % W_{o,f}^{after}
    Of = numel(W_after);
    Wof_a_vec = zeros(Of,1);
    for o = 1:Of
        Wof_a_vec(o) = Wof_fun(W_after{o}, alphaJ_all{f}{o}, params.beta);
    end
    Wof_after_all{f} = Wof_a_vec;
end

% ---------- aggregate \bar{W} ----------
SW = 0.0;
for f = 1:F
    aOf = alphaOf_all{f};
    if isempty(aOf), continue; end
    aF = alphaF(f); Wof = Wof_after_all{f};
    for o = 1:numel(Wof)
        SW = SW + aF * aOf(o) * ( Wof(o)^(1+zeta) );
    end
end
Wbar_imp = SW^(1/(1+zeta));

% ---------- aggregate \bar{Y} ----------
theta = params.theta;
sF = sum( (alphaF(:).^(1/theta)) .* (Yf_after_all(:).^((theta-1)/theta)) );
Ybar_imp = sF^(theta/(theta-1));

% ---------- aggregate \bar{N} via mapping ----------
%{
SN = 0.0;
for f = 1:F
    aOf = alphaOf_all{f};
    if isempty(aOf), continue; end
    aF_inv = alphaF(f)^(-1/zeta);
    Wof    = Wof_after_all{f};
    for o = 1:numel(Wof)
        Nof_a = (Ybar0 * gY) * (zbar0 * S) * ( Wof(o)^zeta );
        SN    = SN + aF_inv * (aOf(o)^(-1/zeta)) * ( Nof_a^((zeta+1)/zeta) );
    end
end
Nbar_imp = SN^(zeta/(zeta+1));
%}
Nbar_imp = Phi^(1/(1+sigma*phi)) * Wbar_imp^(phi*(1-sigma)/(1+sigma*phi)); 
end
