function kint_t = compute_correlation_integral_rate(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;
denom = mean(dh_A.*dh_A);
for t = 1:(tmax+1)
  Cdt(t) = mean(dh_A(1:T-t) .* dh_A(1+t:T)) / denom;
end
Cdt = [1 Cdt];

% Compute correlation integral.
kint_t = zeros(1,tmax);
kint_t(1) = 1.0 / (dt * 0.5 * (Cdt(1) + Cdt(2)));
for t = 2:tmax
  kint_t(t) = 1.0 / ((1.0/kint_t(t-1)) + dt * 0.5 * (Cdt(t) + Cdt(t+1)));
end

return
