function MalvestioPRE2017Example(signal1, signal2, Q, s_div, N_wo, dchoice, k, varargin)
%% MalvestioPRE2017Example v 1.0 02.10.2018
% It shows an example of computation of the measure L with signals from
% Hindmarsh-Rose dynamics and a suited choice of parameter.
% The parameter and signals here are taken from the paper [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.
% For further information, please refer to it, and to MalvestioPRE2017Readme.pdf.
% The execution of this source code can take relatively long
%
% signal1: times of the spikes in the first signal, vector of (1, N1), where N1 is the number of spikes in the first signal.
% signal2: analogous to signal1 but for the second spike train.
% Q: total recording time in the same unit as the spike times.
% s_div: factor of overlap between windows (q/s).
% N_wo: number of windows without overlap. The total number of windows considered is N_w = N_wo*s_div. The computation time will be very long if N_w is too large.
% dchoice: allow to select the spike train distance
% 1 : ISI distance
% 2 : SPIKE distance
% 3 : Adaptive ISI distance
% 4 : Adaptive SPIKE distance
% k: number of nearest neighbors
if nargin < 1
load MalvestioPRE2017Example024
% Parameters as in the paper [1]:
Q = 400000;
s_div = 5;
N_wo = 400;
dchoice = 3;
k = 10;
end
tic
% Define parameters
disp('Defining parameters');
q = floor(Q/N_wo);
s = q/s_div;
W = s_div;
% The value of Wex is set adaptively as the minimum number of windows that we have to neglect in order not to consider windows with empty neighbor windows.
% Minimum discard at the beginning:
W1 = floor(max(signal1(1), signal2(1))/s) +1 ;
% Minimum discard at the end:
W2 = floor((Q-q)/s) - floor((min(signal1(end), signal2(end))-q)/s);
W_ex = max(W1, W2);
% Calculate distance matrices and their corresponding rank matrices HP_1 and HP_2
disp('Calculating distance matrices and the corresponding ranks of individual windows');
HP_1 = MalvestioPRE2017Hpointprocess(signal1,Q,q,s,W,W_ex,dchoice);
HP_2 = MalvestioPRE2017Hpointprocess(signal2,Q,q,s,W,W_ex,dchoice);
% Now all ranks matrices can be treated the same and stacked in array
Hranks(:,:,1) = HP_1;
Hranks(:,:,2) = HP_2;
% Calculation of the actual measure L
disp('Calculation of the actual measure L');
[Lj, L] = MalvestioPRE2017FromHtoL(Hranks,k,W);
disp('In L(i,j) we have L(signal i | signal j)');
disp(L);
toc