function [tau,v_AMI]=AMI_Stergiou(data,L,bins)

% [tau,v_AMI]=AMI(data,L,bins)
% inputs    - data, column oriented time series
%           - L, maximal lag to which AMI will be calculated
%           - bins, number of bins to use in the calculation, if empty an
%             addaptive formula will be used
% outputs   - tau, first minimum in the AMI vs lag plot
%           - v_AMI, vector of AMI values and associated lags
%
% Remarks
% - This code uses average mutual information to find an appropriate lag
%   with which to perform phase space reconstruction. It is based on a
%   histogram method of calculating AMI.
% - In the case a value of atu could not be found before L the code will
%   automatically re-execute with a higher value of L, and will continue to
%   re-execute up to a ceiling value of L.
%
% Future Work
% - None currently.
%
% Mar 2015 - Modified by Ben Senderling, email unonbcf@unomaha.edu
%          - Modified code to output a plot and notify the user if a value
%            of tau could not be found.
% Sep 2015 - Modified by Ben Senderling, email unonbcf@unomaha.edu
%          - Previously the number of bins was hard coded at 128. This
%            created a large amount of error in calculated AMI value and
%            vastly decreased the sensitivity of the calculation to changes
%            in lag. The number of bins was replaced with an adaptive
%            formula well known in statistics. (Scott 1979
%          - The previous plot output was removed.
% Oct 2017 - Modified by Ben Senderling, email unonbcf@unomaha.edu
%          - Added print commands to display progress.
%%

N=length(data);
epsilon = eps;      %or use epsilon = 1e-10;

% set number of bins
if isempty(bins)
    bins=ceil(range(data)/(3.49*std(data)*N^(-1/3))); % Scott 1979
%     fprintf('number of bins = %.0f\n',bins)
end

data = data - min(data); % make all data points positive
data = 1+ floor(data/(max(data)/(bins-epsilon))); % scaling the data

v=zeros(L,1); % preallocate vector
overlap=N-L;

increment= 1/overlap;
one = ones(overlap,1); %create a column vector with all elements being one

% MUTUAL INFORMATION (not used, do not know why, BS ))
% I (time_lag) = sum [ p(x(t), x(t + time_lag))*log[(p(x(t),p(x + time_lag))/p(x(t))*p(x(t+time_lag))]
%find probability p(x(t))= pA

pA = sparse (data(1:overlap),one,increment);  
%e.g. when overalp = N+1-L = 6001+1-32= 5970, max(data(1:overlap))=129, 
%creating a histogram with (129-1) bins
% sum(pA)= 1 --> 100 % in total 

tau=[];
v=zeros(L+1,2);

fprintf('AMI: 00%%')

for lag = 0: L % used to be from 0:L-1 (BS)
    
    v(lag+1,1)=lag;
    
    %find probablity p(x(t+time_lag))=pB, sum(pB)=1
    pB = sparse(one, data(1+lag:overlap+lag), increment);
    %find joing probability p(A,B)=p(x(t),x(t+time_lag))
    pAB = sparse(data(1:overlap),data(1+lag:overlap+lag),increment);
    [A, B, AB]=find(pAB);
    v(lag+1,2)=sum(AB.*log2(AB./(pA(A).*pB(B)')));  %Average Mutual Information
    
    fprintf(repmat('\b',1,4));
    msg=sprintf('%3.0d',floor(lag/L*100));
    fprintf([msg '%%']);
    
end

tau=[];
for i=2:length(v)-1
    if (v(i-1,2)>=v(i,2))&&(v(i,2)<=v(i+1,2))
        tau(end+1,:)=[v(i,1) v(i,2)];
    end
end

v_AMI=v;

% Catch if no value of tau is found
if isempty(tau)
    if L*1.5>length(data/3)
        tau=9999;
    else
        fprintf('max lag needed to be increased for AMI_Stergiou\n')
        [tau,v_AMI]=AMI_Stergiou(data,floor(L*1.5),[]);
    end
end


