function [kim_t, x_ddag] = compute_implied_rate_automatic(x_t, x_ddag_min, x_ddag_max, dt)
% Compute implied rate estimate as a function of observation time.
%
% tmax is automatically selected to be tmax = 2 * (1/k).
% The optimal dividing surface x_ddag is automatically guessed within the specified range.
% The optimal dividing surface x_ddag minimizes the instantaneous crossing rate estimate.
%
% [kim_t, x_ddag] = compute_implied_rate_automatic(x_t, x_ddag_min, x_ddag_max, dt)
%
% ARGUMENTS
%   x_t - timeseries
%   x_ddag_min - minimum value to try for dividing surface separating A and B states
%   x_ddag_max - maximum value to try for dividing surface separating A and B states
%   dt - time interval between observations in x_t
%
% RETURNS
%   kim_t - implied rate as a function of lag time, up to time tmax
%   x_ddag - guess at optimal dividing surface.

% Find optimal dividing surface.
ntry = 100; % number of dividing surfaces to try
x_ddag_try = linspace(x_ddag_min, x_ddag_max, ntry);
k_i = zeros(ntry,1);
for i = 1:ntry
  x_ddag = x_ddag_try(i);
  k_i(i) = compute_crossing_rate(x_t, x_ddag, dt);
end
[k_cross, i] = min(k_i);
x_ddag = x_ddag_try(i);

% 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.
dh_A = h_A - PA;
denom = mean(dh_A.*dh_A);
kim_t = zeros(T,1);

t = 1;
Cdt = mean(dh_A(1:T-t) .* dh_A(1+t:T)) / denom;
kim_t(1) = - (t*dt).^(-1) .* log(Cdt);

while (2*(1.0/kim_t(t)) > t*dt)
  t = t + 1;
  Cdt = mean(dh_A(1:T-t) .* dh_A(1+t:T)) / denom;
  kim_t(t) = - (t*dt).^(-1) .* log(Cdt);
end

% Truncate.
kim_t = kim_t(1:t);

return
