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 & spiketimes0+wc1*s & spiketimes