%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% AMI wrapper algorithm
function [tau] = modified_estimateTauAMI(numbins, omega_uniform, tauPlot)

 omega_bins    = linspace(min(omega_uniform), max(omega_uniform), numbins)';

 % initialize AMI minimum search; then calculate AMIs until minimum found (greedy)
 av_I(1) = find_av_I(omega_uniform,1, omega_bins, numbins);
 av_I(2) = find_av_I(omega_uniform,2, omega_bins, numbins);
 t       = 2;
 while (av_I(t-1) > av_I(t))
    t       = t+1;
    av_I(t) = find_av_I(omega_uniform, t,omega_bins, numbins);
 end;
%  for t=1:50
%  %   t       = t+1;
%     av_I(t) = find_av_I(omega_uniform, t,omega_bins, numbins);
%  end;
 % stopped loop when first minimum found
 if (tauPlot > 0) figure(10); plot([1:t],av_I,'o'); end;
 tau = t-1;
 %tau=12;

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% algorithm for finding average mutual information (AMI) of vector
function av_I=find_av_I(omega_uniform,T_lag,omega_bins,numbins)

 % individual probabilities
 omega_delay = [omega_uniform(1:end-T_lag), omega_uniform(T_lag+1:end)]; 
 omega_hist  = hist(omega_delay,omega_bins); 
 omega_prob  = omega_hist/sum(omega_hist(:,1),1);

 % tabulate joint distribution
 binomega   = 1+floor((numbins-1).*(omega_uniform-min(omega_uniform))/(max(omega_uniform)-min(omega_uniform)));
 joint_dist = zeros(numbins);                  % initialize empty joint distribution matrix
 for ii=1:size(binomega,1)-T_lag               % count co-occurrences
     joint_dist(binomega(ii),binomega(ii+T_lag)) = joint_dist(binomega(ii),binomega(ii+T_lag)) +1;
 end
 joint_dist=joint_dist/(sum(sum(joint_dist))); % normalize joint distribution

 % calculate the AMI
 warning off
 prob_ratio = joint_dist./(omega_prob(:,1)*omega_prob(:,2)');
 warning on
 prob_ratio(find(joint_dist==0)) = 1;
 prob_ratio(isinf(prob_ratio))   = 1;
 mutual_information              = log2(prob_ratio);

 % return result
 av_I=sum(sum(joint_dist.*mutual_information));
