function [L] = LMultiBig(Datarec,m,tau,krec,W)
%
% Calculates Interdependence measures
% H, S:   J. Arnhold, K. Lehnertz, P. Grassberger, and C. E. Elger,
%         Physica D 134, 419 1999.
% N:      R. Quian Quiroga, A. Kraskov, T. Kreuz, and P. Grassberger,
%         Phys. Rev. E 65, 041903 2002.
% M:      defined independently in (A) R. G. Andrzejak, A. Kraskov, H. St�gbauer, F. Mormann, and
%                                      T. Kreuz, Phys. Rev. E 68, 066202 2003.
%                                  (B) H. Kantz and T. Schreiber, Nonlinear Time Series Analysis
%                                      Cambridge University Press, Cambridge, U.K., Second edition 2003, Chapter 14
% L:     D. Chicharro, A. Ledberg, R.G. Andrzejak in preparation
%
% The data is  expected in Datarec with dimension (N*,2). N* is the length
% of the original time series. The first component should contain X, the
% second should contain Y. m is the embedding dimension and tau is the time
% delay used for the delay reoncstruction. krec is an array of numbers of
% nearest neighbors k. The algorithm can return results for a whole range
% of k in just one run. W is the Theiler correction. W temporally nearest
% neighbors will be excluded from the distance matrix. For W = 0, only the
% diagonal will be excluded.
% The results will be found in H organized in an (5,2,length(krec)) array.
% The first axis is for [H, L, S, M, N]
% The second axis is for the directions IMPORTANT it returns [(Y|X) (X|Y)]
% The third axis is for the different numbers of nearest neighbors [krec]
%
%
%   The original code has been optimized to speed up the calculations for a
%   large number of pairs of signals as we have in PRE, 99.1 (2019) 012319.
%  
%


% Determine the number of data points and substract the embedding window
PHrec = [];
NM = size(Datarec);
if ~isempty(PHrec)
    N = NM(1)-max(PHrec);
else
    N = NM(1);
end
D = NM(2);
embwin = (m-1)*tau;
Neff = N -embwin;

% Ncorrected(counter) is the effective number of distances in row/colum
% #counter of the distance matrix. Here the Theiler correction that skipps
% W data points has to be taken into account

Ncorrected = zeros(1,Neff);
Ncorrected(1:Neff) = Neff-(W*2+1);
for counter = 1:W
    Ncorrected(counter) = Neff-W-(counter);
    Ncorrected(Neff-counter+1) = Neff-W-(counter);
end


% Allocate matrices
Hdistances =zeros(Neff,Neff,'single');
Hranks = zeros(Neff,Neff,D,'int16');
Hindexes = zeros(Neff,Neff,D,'int16');

% Determine distance matrices
for counter3 = 1:D
    for counter1 = embwin+1:embwin+1+tau
        for counter2 = counter1+1+W:N
            distA = 0;
            
            for index = -embwin:tau:0;
                distA = distA+(Datarec(counter1+index,counter3)-Datarec(counter2+index,counter3))^2; % X
            end
            Hdistances(counter1-embwin,counter2-embwin)=distA;
            Hdistances(counter2-embwin,counter1-embwin)=distA;
            
        end
        
    end
    
    for counter1 = embwin+2+tau:N
        for counter2 = counter1+1+W:N
            Hdistances(counter1-embwin,counter2-embwin)= Hdistances(counter1-embwin-tau,counter2-embwin-tau)-(Datarec(counter1-embwin-tau,counter3)-Datarec(counter2-embwin-tau,counter3))^2+(Datarec(counter1,counter3)-Datarec(counter2,counter3))^2;
            Hdistances(counter2-embwin,counter1-embwin)= Hdistances(counter1-embwin,counter2-embwin);
        end
    end
    
    %    Hsums(:,counter3) = sum(Hdistances(:,:,counter3))./Ncorrected;
    highnumber = 1e99;
    Hdistances(Hdistances==0)=highnumber;
    [~, Hindexes(:,:,counter3)] = sort(Hdistances(:,:));
    N1 = 1:Neff;
    for s = 1:Neff
        Hranks(Hindexes(:,s,counter3),s,counter3) = N1;
        
    end
    
    
end
% whos Hsums
% overwrite zeros that were skipped due to the Theiler correction with an
% arbitrary high number. IMPORTANT: This number must be higher than the
% highest distance expected in the distance matrix
% highnumber = 1e99;
% Hdistances(Hdistances==0)=highnumber;
% 
% for counter3 =1:D
%     % Calculate the matrices of ranks
%     % Note that were here use the same way to determine the ranks like Matlab
%     % does for example in the function ranksum.m
%     
%     
%     [~, Hindexes(:,:,counter3)] = sort(Hdistances(:,:,counter3));
%     N1 = 1:Neff;
%     for s = 1:Neff
%         Hranks(Hindexes(:,s,counter3),s,counter3) = N1;
%         
%     end
% end

%Hj= zeros(D,D,length(krec));
%Sj= zeros(D,D,length(krec));
Lj= zeros(D,D,length(krec));
%Mj= zeros(D,D,length(krec));
%Zj= zeros(D,D,length(krec));

help1 = zeros(1,Neff);
for counter = 1:Neff
    help1(counter) = (Ncorrected(counter)+1)/2;
end

kc = 0;
for k = krec
    help2 = (k+1)/2;
    kc = kc +1;
    for j = 1:Neff
        for counter3 = 1:D
            %Rk = sum(Hdistances(Hindexes(1:k,j,counter3),j,counter3))/k;
            for counter4 = [1:counter3-1 counter3+1:D]
                %Rkcond = sum(Hdistances(Hindexes(1:k,j,counter4),j,counter3))/k;
                %Hj(counter3,counter4,kc)=  Hj(counter3,counter4,kc)+log(Hsums(j,counter3)/Rkcond);
                Lj(counter3,counter4,kc) = Lj(counter3,counter4,kc)+(help1(j)-sum(Hranks(Hindexes(1:k,j,counter4),j,counter3))/k)/(help1(j)-help2);
                %Sj(counter3,counter4,kc) = Sj(counter3,counter4,kc)+Rk/Rkcond;
                %Mj(counter3,counter4,kc) = Mj(counter3,counter4,kc)+(Hsums(j,counter3)-Rkcond)/(Hsums(j,counter3)-Rk);
                %Zj(counter3,counter4,kc) = Zj(counter3,counter4,kc)+(Hsums(j,counter3)-Rkcond)/(Hsums(j,counter3));
            end
        end
    end
end

% Finally, take the mean across all reference points

L = Lj/Neff;


end