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;