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;
