function kRF_t = compute_reactive_flux(x_t, x_ddag, dt, tmax)
% Compute reactive flux as a function of observation time.
%
% kRF_t = compute_reactive_flux(x_t, x_ddag, tmax)
%
% ARGUMENTS
%   x_t - timeseries
%   x_ddag - dividing surface separating A and B states
%   dt - time interval between observations in x_t
%   tmax - maximum lag time to compute implied rate
%
% RETURNS
%   kRF_t - reactive flux as a function of lag time, up to time tmax

% Determine length of timeseries.
T = length(x_t);

% Quantize using dividing surface to define h_A and h_B functions.
h_A = (x_t <= x_ddag);
h_B = (x_t > x_ddag);

% Compute stationary probabilities.
PA = mean(h_A);
PB = mean(h_B);

% Compute normalized time-correlation function for dh_A.
Cdt = zeros(1,tmax+1);
dh_A = h_A - PA;
for t = 1:(tmax+1)
  Cdt(t) = mean(dh_A(1:T-t) .* dh_A(1+t:T));
end
Cdt = [1 Cdt];

% Compute reactive flux correlation function
kRF_t = - (Cdt(2:tmax+1) - Cdt(1:tmax)) / dt;
kRF_t = kRF_t / mean(dh_A.*dh_A); % assuming <h_A(0) h_A(t)> ~ <h_A(0) h_A(0)>
%kRF_t = kRF_t ./ Cdt(1:tmax); % avoiding this assumption

return
