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 