function L = AndrzejakKreuzL(DArray,k,W)

% For further information, please refer to AndrzejakKreuzReadme.pdf and to
% Andrzejak RG, Kreuz T: Characterizing unidirectional couplings between point processes and flows. EPL, 96 (2011) 50012.


% Determine the number of data points and substract the embedding window
N = size(DArray,1);
Neff = N;
D = size(DArray,3);

% 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,D,'single');
Hranks = zeros(Neff,Neff,D,'int16');
Hindexes = zeros(Neff,Neff,D,'int16');

% Determine distance matrices
for counter3 = 1:D
    Hdistances(:,:,counter3)=DArray(:,:,counter3);          
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


    [dummy, Hindexes(:,:,counter3)] = sort(Hdistances(:,:,counter3));
    % dummy is not used. Depending on your matlab version you replace it by
    % the tilde symbol.
    N1 = 1:Neff;
    for s = 1:Neff
        Hranks(Hindexes(:,s,counter3),s,counter3) = N1;

    end
end

Lj= zeros(D,D);

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



help2 = (k+1)/2;
for j = 1:Neff
    for counter3 = 1:D
        for counter4 = [1:counter3-1 counter3+1:D]
            Lj(counter3,counter4) = Lj(counter3,counter4)+(help1(j)-sum(Hranks(Hindexes(1:k,j,counter4),j,counter3))/k)/(help1(j)-help2);
        end
    end
end


% Finally, take the mean across all reference points

L = Lj/Neff;

