
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%          Motif Entropy
%                      A program for calculating the entropy
%                      of motifs in a time series.
%                      Developed by M. Zanin.
%
%                      V1.0
%
%                      Input arguments:
%                          TS: (t x s) time series, composed of t samples
%				of s different time series
%                          windowSize: window size for calculating correlations
%                          synchThr: correlation threshold [0, 1]
%
%                      Outputs:
%                          normMotifEntropy: normalized motif entropy
%
%                      Please cite the following paper:
%                          M. Zanin, D. Papo,
%                          Characterizing Motif Dynamics of Electric Brain Activity Using Symbolic Analysis
%                          Entropy 2014, 16(11), 5654-5667
%


function [normMotifEntropy] = MotifEntropy(TS, windowSize, synchThr)

numSamples = size(TS, 1);
numWindows = floor(numSamples ./ windowSize);

motifEvolution = zeros( numWindows, 1 );

for offset = 1 : numWindows
    
    off1 = (offset - 1) * windowSize + 1;
    off2 = (offset) * windowSize;

    subTS = TS( off1 : off2, : );
    
    corrMatrix = zeros(3, 3);
    
    for n1 = 1 : 3
        for n2 = 1 : 3
            
            if n1 == n2
                continue;
            end
            
            corrMatrix(n1, n2) = corr( subTS(:, n1), subTS(:, n2) );
            
        end
    end
            
    corrMatrix = abs( corrMatrix ) > synchThr;
    
    motifNumber = bi2de( corrMatrix(:)' );
    motifEvolution( offset ) = motifNumber;

end

motifEvolution( motifEvolution == 0 ) = 1;
motifEvolution( motifEvolution == 10 ) = 2;
motifEvolution( motifEvolution == 68 ) = 3;
motifEvolution( motifEvolution == 78 ) = 4;
motifEvolution( motifEvolution == 160 ) = 5;
motifEvolution( motifEvolution == 170 ) = 6;
motifEvolution( motifEvolution == 228 ) = 7;
motifEvolution( motifEvolution == 238 ) = 8;

motifProb = zeros(5, 1);
for k = 1 : 8
    motifProb(k) = sum( motifEvolution == k ) ./ numWindows;
end

motifEntropy = 0;
for k = 1 : 8
    if motifProb(k) == 0
        continue;
    end
    
    motifEntropy = motifEntropy - motifProb(k) * log2( motifProb(k) );
end

normMotifEntropy = motifEntropy ./ log2(8);

end





