function tau = integrated_autocorrelation_time(A_n)
% Compute the integrated autocorrelation time of a timeseries.
%
%  REQUIRED ARGUMENTS  
%    A_n (numpy array) - A_n[n] is nth value of timeseries A.  Length is deduced from vector.
%
%  OPTIONAL ARGUMENTS
%    B_n (numpy array) - B_n[n] is nth value of timeseries B.  Length is deduced from vector.
%       If supplied, the cross-correlation of timeseries A and B will be estimated instead of the
%       autocorrelation of timeseries A.  
%    fast (boolean) - if True, will use faster (but less accurate) method to estimate correlation
%       time, described in Ref. [1] (default: False)
%    mintime (int) - minimum amount of correlation function to compute (default: 3)
%       The algorithm terminates after computing the correlation time out to mintime when the
%       correlation function furst goes negative.  Note that this time may need to be increased
%       if there is a strong initial negative peak in the correlation function.
%
%  RETURNS
%    g is the estimated statistical inefficiency (equal to 1 + 2 tau, where tau is the correlation time).
%       We enforce g >= 1.0.
%
%  NOTES 
%    The same timeseries can be used for both A_n and B_n to get the autocorrelation statistical inefficiency.
%    The fast method described in Ref [1] is used to compute g.
%
%  REFERENCES  
%    [1] J. D. Chodera, W. C. Swope, J. W. Pitera, C. Seok, and K. A. Dill. Use of the weighted
%    histogram analysis method for the analysis of simulated and parallel tempering simulations.
%    JCTC 3(1):26-41, 2007.

% Set defaults.
fast = 1;
mintime = 3;

% Get the length of the timeseries.
N = length(A_n);

% Initialize autocorrelation estimate with uncorrelated value.
tau = 0.0;
    
% Compute mean of timeseries.
mu_A = mean(A_n);

% Make temporary copies of fluctuation from mean.
dA_n = A_n - mu_A;

% Compute estimator of covariance of (A,B) using estimator that will ensure C(0) = 1.
sigma2_AA = mean(dA_n .* dA_n); % standard estimator to ensure C(0) = 1

% Accumulate the integrated correlation time by computing the normalized correlation time at
% increasing values of t.  Stop accumulating if the correlation function goes negative, since
% this is unlikely to occur unless the correlation function has decayed to the point where it
% is dominated by noise and indistinguishable from zero.
t = 1;
increment = 1;
while (t < N-1)
    % compute normalized fluctuation correlation function at time t
    C = mean(dA_n(1:(N-t)) .* dA_n((1+t):N)) / sigma2_AA;
    
    % Terminate if the correlation function has crossed zero and we've computed the correlation
    % function at least out to 'mintime'.
    if ((C <= 0.0) && (t > mintime))
      break
    end
    
    % Accumulate contribution to the statistical inefficiency.
    tau = tau + C * (1.0 - t/N) * increment;

    % Increment t and the amount by which we increment t.
    t = t + increment;

    % Increase the interval if "fast mode" is on.
    if fast
        increment = increment + 1;
    end
end

% tau must be nonnegative
if (tau < 0.0)
    tau = 0.0;
end
   
% Return the computed statistical inefficiency.
return


