function [Lj, L] = MalvestioPRE2017FromHtoL(Hranks,k,W)
%% MalvestioPRE2017FromHtoL v 2.0 18.10.2018
% Computes the values of L between the all the pairs of signals.
% Lj are the values in the right side of Eq. 2 of the paper [1], before the average.
%
% Hranks: for each signal w, Hranks(i,j,w) is g_{i,j} for the signal w.
% k: number of neighbours considered
% W: windows excluded due to Theiler correction
%
% For further information, please refer to
% MalvestioPRE2017Readme.pdf
% and to [1]:
% Malvestio I, Kreuz T, Andrzejak RG:
% Robustness and versatility of a nonlinear interdependence method for directional coupling detection from spike trains.
% PRE, 96.2 (2017) 022203.
N_node = size(Hranks,3); % number of nodes between which we are computed L
% From the matrix with the rank of the distance, to the matrix with the indexes of the corresponding window:
[~, Hindexes] = sort(Hranks);
% Determine the number of windows of the matrixes
Neff = size(Hranks,1);
% Compute the values for the normalization of L:
help1 = (Neff-(W*2+1))/2 * ones(1,Neff);
% Correct this value taking into account of the effect of Theiler correction
for counter = 1:W
help1(counter) = (Neff-W-(counter)+1)/2;
help1(Neff-counter+1) = (Neff-W-(counter)+1)/2;
end
% Compute another value for the normalization of L:
help2 = (k+1)/2;
L = zeros(N_node, N_node);
Lj = zeros(N_node, N_node, Neff);
% Computation of Lj from signal node1 to signal node2 and vice versa
for node1 = 1: N_node
for node2 = (node1+1): N_node %%%starting from node1 +1 we avoid to compute between a node and itself, and between pairs already computed
for j = 1:Neff
Lj(node1,node2, j) = (help1(j)-sum(Hranks(Hindexes(1:k,j,node2),j,node1))/k)/(help1(j)-help2);
Lj(node2,node1, j) = (help1(j)-sum(Hranks(Hindexes(1:k,j,node1),j,node2))/k)/(help1(j)-help2);
end
end
end
% Finally, L is the mean of Lj across all reference points
L(:,:) = sum(Lj,3) /Neff;