function s = fastSampen(y,m,r)
% INPUTS
%  y     num_var x num_samples
%  m     template length
%  r     radius
% OUTPUTS
%  s     Sample Entropy
%
% Written by Shamim Nemati <shamim.nemati@alum.mit.edu>
%
% 01-31-2018 : modified by Giulia Da Poian <giulia.dap@gmail.com>, see
%              comments in the code 
%              NOTE : Sample entropy quantifies the likelihood that a 
%              sequence of m consecutive data points that matches another
%              sequence of the same length (match within a tolerance of r) 
%              will still match the other sequence when their length is 
%              increased of one sample (sequences of length m + 1); 
%              References : 
%              1. Richman JS, Moorman JR. Physiological time-series 
%              analysis using approximate entropy and sample entropy. 
%              American Journal of Physiology-Heart and Circulatory 
%              Physiology. 2000 Jun 1;278(6):H2039-49.
%              2. Humeau-Heurtier A. The multiscale entropy 
%              algorithm and its variants: A review. Entropy. 2015 May 
%              12;17(5):3110-23.
% LICENSE:    
%       This software is offered freely and without warranty under 
%       the GNU (v3 or later) public license. See license file for
%       more information

% Ensure y is a row vector
if size(y, 1) > size(y,2)
    y = y';
end

xx = convert_to_lagged_form(y, m)';   % Giulia: replaced m-1 with m to match SampEn definition
Dxx = pdist(xx,'chebychev');

yy = convert_to_lagged_form(y, m+1)'; % Giulia: replaced m with m+1 to match SampEn definition
Dyy = pdist(yy,'chebychev');

A = mean( Dxx < r ) ;
B = mean( Dyy < r );

s = -log(B/A);


function yy = convert_to_lagged_form(y, k)
% Create an observation vector yy(:,t) containing the last k values of y, newest first
% e.g., k=2, y = (a1 a2 a3)     yy  = a2 a3
%                (b1 b2 b3)           b2 b2
%                                     a1 a2
%                                     b1 b2
[s, T] = size(y);
bs = s*ones(1,k);
yy = zeros(k*s, T-k+1);
for i=1:k,     yy(block(i,bs), :) = y(:, k-i+1:end-i+1); end

function sub = block(blocks, block_sizes)
% BLOCK Return a vector of subscripts corresponding to the specified blocks.
% sub = block(blocks, block_sizes)
%
% e.g., block([2 5], [2 1 2 1 2]) = [3 7 8].
blocks = blocks(:)';
block_sizes = block_sizes(:)';
skip = [0 cumsum(block_sizes)];
start = skip(blocks)+1;
fin = start + block_sizes(blocks) - 1;
sub = [];
for j=1:length(blocks)
    sub = [sub start(j):fin(j)];
end