function HP = MalvestioPRE2017Hpointprocess(spiketimes,Q,q,s,W,W_ex,dchoice)

%% MalvestioPRE2017Hpointprocess  v 1.0 02.10.2018
% Adaptation of AndrzejakKreuzDpointprocess with the implemention of spike
% train distance made by Eero Satuvuori (http://wwwold.fi.isc.cnr.it/users/thomas.kreuz/Source-Code/cSPIKE.html ).
% It calculates the distances between every pairs of windows of the same spike train, and from this the corresponding
% matrix of the ranks.
% In particular, HP(i,j) = g_{i,j}, introduced in the paper [1] before Eq. 1. 
%
% spiketimes: 1D vector containing times of spikes
% Q: interval of time during which two dynamics are simultaneously measured.
% q: breakdown of Q into smaller overlaping intervals (q < Q).
% s: step size of the overlaping intervals q  (s < q).  
% s and q need to be choose in order to have (Q-q)/s as an integer(otherwise the very end of the signal may be ignored).
% W: Theiler Window, to avoid comparison between neighboring windows
% W_ex: number of windows excluded at the very beginning and at the very end.
% dchoice allow to select the distance: 
%   1 : ISI distance 
%   2 : SPIKE distance
%   3 : Adaptive ISI distance
%   4 : Adaptive SPIKE distance
%
% 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.

InitializecSPIKE; 
N=fix((Q-q)/(s)+1);
DP=zeros(N,N);

if dchoice == 3 || 4 
    duration_total = Q; 
    spikes_total1 = spiketimes(spiketimes>0 & spiketimes< duration_total); 
    spikemat_total{1} = spikes_total1;
    spikemat_total{2} = spikes_total1;
    STS_total = SpikeTrainSet(spikemat_total, 0, duration_total);
    [threshold, ~] = giveTHR(STS_total);
end

for wc1=W+1:N-1   % calculate distance profile
    winspikes1=spiketimes(spiketimes>0 & spiketimes<Q-wc1*s);
    winspikes2=spiketimes(spiketimes>0+wc1*s & spiketimes<Q)-wc1*s;
    duration=Q-wc1*s;

    spikemat{1} = winspikes1;
    spikemat{2} = winspikes2;

    STS = SpikeTrainSet(spikemat, 0, duration);
  
    for wc2=1:N-wc1    % calculate the actual distance value that is the average into the windows
        wmin=(wc2-1)*s;
        wmax=wmin+q;
        switch  dchoice
            case {1}  % ISI distance
                DP(wc1+wc2,wc2) = STS.ISIdistance(wmin, wmax);  
            case {2} % SPIKE distance 
                DP(wc1+wc2,wc2) = STS.SPIKEdistance(wmin, wmax);  
            case {3} % adaptative ISI distance
                 DP(wc1+wc2,wc2) = STS.AdaptiveISIdistance(wmin, wmax, threshold);   %% we need to use the same threshold for all the windows
            case {4 } % adaptive SPIKE distance
                 DP(wc1+wc2,wc2) = STS.AdaptiveSPIKEdistance(wmin, wmax, threshold); %% we need to use the same threshold for all the windows 
        end
        DP(wc2,wc1+wc2) = DP(wc1+wc2,wc2);     %makes the matrix simmetric      
    end 
end

N = size(DP,1); %total number of windows 

% Discard W_ex windows at the beginning and at the end (in order to avoid edge effect)

begg = W_ex +1;
endd = N - W_ex;
DP = DP(begg:endd, begg:endd);
Neff = size(DP,1); 

% Calculation of the corresponding rank:
Hdistances = zeros(Neff,Neff, 'single');
HP = zeros(Neff,Neff, 'int16'); 
Hindexes = zeros(Neff,Neff, 'int16'); 

Hdistances = DP; 
clear DP;
% 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;

[~, Hindexes] = sort(Hdistances);
clear Hdistances;
[~, HP] = sort(Hindexes); 