function H = HSLMNCom(Datarec,m,tau,krec,W)
%
% Calculates the following 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, T. Kreuz, Phys. Rev. E 68, 066202, 2003.
% and in
% H. Kantz and T. Schreiber, Nonlinear Time Series Analysis Cambridge University Press, Cambridge, U.K., Second edition 2003, Chapter 14
% D. Chicharro and R. G. Andrzejak Reliable detection of directional couplings using rank statistics. Physical Review E, 80, 026217, 2009
%
%
%
% The data must be given 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.
%
% The parameter m is the embedding dimension and tau is the time
% delay used for the delay reoncstruction.
%
% The parameter 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.
%
% The parameter 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 are output in H organized in an (5,2,length(krec)) array.
% The first axis comprises the different measures [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]
% Algorithm Start
% Determine the number of data points and substract the embedding window.
N = size(Datarec);
N = N(1);
embwin = (m-1)*tau;
Neff = N -embwin;
% Allocate matrices.
Hdistances1 =zeros(Neff,Neff);
Hdistances2 =zeros(Neff,Neff);
Hranks1 = zeros(Neff,Neff);
Hranks2 = zeros(Neff,Neff);
% Determine distance matrices.
for counter1 = embwin+1:embwin+tau
for counter2 = counter1+1+W:N
distA = 0;
distB = 0;
for index = -embwin:tau:0;
distA = distA+(Datarec(counter1+index,1)-Datarec(counter2+index,1))^2; % X
distB = distB+(Datarec(counter1+index,2)-Datarec(counter2+index,2))^2; % Y
end
Hdistances1(counter1-embwin,counter2-embwin)=distA;
Hdistances2(counter1-embwin,counter2-embwin)=distB;
end
end
for counter1 = embwin+1+tau:N
for counter2 = counter1+1+W:N
Hdistances1(counter1-embwin,counter2-embwin)= Hdistances1(counter1-embwin-tau,counter2-embwin-tau)-(Datarec(counter1-embwin-tau,1)-Datarec(counter2-embwin-tau,1))^2+(Datarec(counter1,1)-Datarec(counter2,1))^2;
Hdistances2(counter1-embwin,counter2-embwin)= Hdistances2(counter1-embwin-tau,counter2-embwin-tau)-(Datarec(counter1-embwin-tau,2)-Datarec(counter2-embwin-tau,2))^2+(Datarec(counter1,2)-Datarec(counter2,2))^2;
end
end
% Symmetrize the distance matrices.
Hdistances1 = Hdistances1+Hdistances1';
Hdistances2 = Hdistances2+Hdistances2';
% Calculate effective number of distances
% Ncorrected(counter) is the number of distances contained in the distance matrix' row/colum
% with index counter. Here, the Theiler correction W 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;
Hsums1 = sum(Hdistances1)./Ncorrected;
Hsums2 = sum(Hdistances2)./Ncorrected;
% 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;
Hdistances1(Hdistances1==0)=highnumber;
Hdistances2(Hdistances2==0)=highnumber;
% Calculate the matrices of ranks.
% Note that were here use the same way to determine the ranks as Matlab
% does for example in the function ranksum.m
[temparr1 Hindexes1] = sort(Hdistances1);
[temparr1 Hindexes2] = sort(Hdistances2);
N1 = 1:Neff;
for s = 1:Neff
Hranks1(Hindexes1(:,s),s) = N1;
Hranks2(Hindexes2(:,s),s) = N1;
end
Hj= zeros(5,2,Neff,length(krec));
% Calculate the actual measures.
kc = 0;
for k = krec
kc = kc +1;
for j = 1:Neff
RkY = sum(Hdistances2(Hindexes2(1:k,j),j))/k;
RkX = sum(Hdistances1(Hindexes1(1:k,j),j))/k;
RcondYX = sum(Hdistances2(Hindexes1(1:k,j),j))/k;
RcondXY = sum(Hdistances1(Hindexes2(1:k,j),j))/k;
Hj(1,1,j,kc)= log(Hsums2(j)/RcondYX);
Hj(1,2,j,kc)= log(Hsums1(j)/RcondXY);
%RNY/Rcond(Y|X)=>H(Y|X)
%RNX/Rcond(X|Y)=>H(X|Y)
Hj(2,1,j,kc) = ((Ncorrected(j)+1)/2-sum(Hranks2(Hindexes1(1:k,j),j))/k)/((Ncorrected(j)+1)/2-(k+1)/2);
Hj(2,2,j,kc) = ((Ncorrected(j)+1)/2-sum(Hranks1(Hindexes2(1:k,j),j))/k)/((Ncorrected(j)+1)/2-(k+1)/2);
% (GN-G(Y|X))/(GN-Gk) => L(Y|X)
% (GN-G(X|Y))/(GN-Gk) => L(X|Y)
Hj(3,1,j,kc) = RkY/RcondYX;
Hj(3,2,j,kc) = RkX/RcondXY;
%RkY/Rcond(Y|X)=>S(Y|X)
%Rk1/Rcond(X|Y)=>S(X|Y)
Hj(4,1,j,kc) = (Hsums2(j)-RcondYX)/(Hsums2(j)-RkY);
Hj(4,2,j,kc) = (Hsums1(j)-RcondXY)/(Hsums1(j)-RkX);
%(RNY-Rcond(Y|X))/(RNY-RkY)=>M(Y|X)
%(RNX-Rcond(X|Y))/(RNX-RkX)=>M(X|Y)
Hj(5,1,j,kc) = (Hsums2(j)-RcondYX)/(Hsums2(j));
Hj(5,2,j,kc) = (Hsums1(j)-RcondXY)/(Hsums1(j));
%(RNY-Rcond(Y|X))/(RNY)=>N(Y|X)
%(RNX-Rcond(X|Y))/(RNX)=>N(X|Y)
end
end
% Finally, take the mean across all reference points
H = squeeze(mean(Hj,3));