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