function k = compute_crossing_rate(x_t, x_ddag, dt)
% Compute rate estimate from number of crossings.
%
% k = compute_crossing_rate(x_t, x_ddag, dt)
%
% 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
%   k - crossing rate


% Get trajectory length.
T = length(x_t);

% Compute indicator functions for A and B.
h_A = (x_t <= x_ddag);
h_B = (x_t > x_ddag);

% Count number of crossings.
h_AB = (h_A(1:end-1) & h_B(2:end)); % A -> B transitions
h_BA = (h_B(1:end-1) & h_A(2:end));
ncross = sum( h_AB | h_BA );

% Estimate rate.
k = double(ncross) / (T * dt);

return
