function res = pesct(dat,wlen,tau,taurange,pval,ensemble,silent)
% PESCT Permutation Ensemble Symbolic Coupling Traces
% The pesct function calculates the symbolic coupling traces using a
% permutation approach. To allow for an analysis of transient events an 
% ensemble of time series can be used to significantly shorten the required
% window size for a time-resolved analysis.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% RES = PESCT(DAT,WLEN,TAU,TAURANGE,PVAL,ENSEMBLE,SILENT) computes the 
% difference RES between the parameters 'T' and 'Tbar' using the 'Symbolic 
% Coupling Traces' method for a data array DAT using the word length WLEN,
% the time lag between consecutive points in the series TAU, and the scalar 
% vector TAURANGE to determine the time lags for the coupling analysis. A
% significance test with significance level PVAL is performed. The
% parameters ENSEMBLE and SILENT determine, whether an ensemble approach is
% used and whether there will be a graphical representation of the results,
% respectively. Depending on the value of ENSEMBLE the result RES will 
% either be a vector (ENSEMBLE==false) or a matrix (ENSEMBLE==true) 
% containing the corresponding values for each time lag specified in 
% TAURANGE (in the second case for each time point).
%
%   'dat'      - (array) [2,N,M] data to analyse, '2' variables (x,y),
%                'N' observations, 'M' realizations (ensemble)
%   'wlen'     - (integer) the desired word length
%   'tau'      - (1x2 array) time steps between letters of the words
%                 for each variable
%   'taurange' - (vector) all time lags to be considered (negative time
%                lags -|tau| represent a shift so that x(t-|tau|)
%                corresponds to y(t))
%   'pval'     - (double) significance threshold
%   'ensemble' - (boolean) use ensemble approach (true) or not (false)
%   'silent'   - (boolean) plot (false) or no plot (true)
%
%
% References:
%
% Wessel N, Suhrbier A, Riedl M, Marwan N, Malberg H, et al. (2009) 
% Detection of time-delayed interactions in biosignals using symbolic 
% coupling traces. Europhys Lett 87: 10004.
%
% Mueller A, Riedl M, Penzel T, Bonnemeier H, Kurths J, et al. (2013) 
% Coupling analysis of transient cardiovascular dynamics. Biomed Tech 
% (Berl) 58: 131–139.
%
% Mueller A, Riedl M, Penzel T, Kurths J, Wessel N (2014) 
% Kardiorespiratorische Koordination und Ensemble-Kopplungsspuren zur 
% ereignisbasierten Charakterisierung kardiovaskulärer Interaktionen 
% während des Schlafes. Somnologie - Schlafforsch und Schlafmedizin 18: 
% 243–251.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% example: bivariate AR-model
% 
% len = 500;
% N = 100;
% M = 200;
% dat = zeros(2,N,M);
% noise = randn(2,len,M);
% for mm = 1:M
%     x = zeros(2,len);
%     x(:,1:2) = randn(2,2);
%     for ii = 3:len
%         x(1,ii) = 0.3*x(1,ii-1) + 0.7*x(2,ii-1) + 0.01*noise(1,ii,mm);
%         x(2,ii) = 0.3*x(2,ii-1) - 0.7*x(1,ii-2) + 0.01*noise(2,ii,mm);
%     end
%     dat(:,:,mm) = x(:,end-N+1:end);
% end
% 
% res = pesct(dat,4,1,-10:10,0.05,1,0);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

d = size(dat);
% variables
vars = d(1);
if vars ~= 2
    error('Please use exactly two variables!');
end
% observations
N = d(2);
% ensemble size
if ensemble
    M = d(3);
else
    M = 1;
end

word_n = factorial(wlen);

wseq = zeros(2,N-(wlen-1)*tau,M);
shift_mat_ind = reshape(0:tau:(wlen-1)*tau,[],1) * ones(1,N-(wlen-1)*tau) +...
    ones(wlen, 1) * reshape(1:(N-(wlen-1)*tau),1,[]);

for ii=1:M
    for jj=1:2
        x = squeeze(dat(jj,:,ii));
        shift_mat = x(shift_mat_ind);
        % sort matrix along rows to build rank orders, equal values retain
        % order of appearance
        [~, sort_ind_mat] = sort(shift_mat,1);
        ind_mat = zeros(size(sort_ind_mat));
        for kk=1:size(ind_mat,2)
            ind_mat(sort_ind_mat(:,kk),kk) = 1:wlen;
        end
        % assign unique number to each pattern (base-wlen number system)
        wseq(jj,:,ii) = wlen.^(0:(wlen-1)) * (ind_mat-1);
    end
end

wseq = squeeze(wseq);

if ensemble
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % ensemble sct
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    T = zeros(N-(wlen-1)*tau,length(taurange)) + NaN;
    Tbar = zeros(N-(wlen-1)*tau,length(taurange)) + NaN;
    
    jj = 0;
    for tt=taurange
        %disp(tt);
        jj = jj + 1;
        % adjust indices according to current time lag
        startv = max(1,1-tt);
        endv = min(N-(wlen-1)*tau,N-(wlen-1)*tau-tt);
        for ii=startv:endv
            % averaging over the ensemble
            % same words --> difference '0'
            % diametral words --> sum 'wlen^wlen-1'
            T(ii,jj) = sum((wseq(1,ii+tt,:)-wseq(2,ii,:)) == 0)/M;
            Tbar(ii,jj) = sum((wseq(1,ii+tt,:)+wseq(2,ii,:)) == wlen^wlen-1)/M;
        end
    end
else
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % standard sct
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    T = zeros(1,length(taurange)) + NaN;
    Tbar = zeros(1,length(taurange)) + NaN;
    jj = 0;
    for tt=taurange
        jj = jj + 1;
        % adjust indices according to current time lag
        startv = max(1,1-tt);
        endv = min(N-(wlen-1)*tau,N-(wlen-1)*tau-tt);
        % averaging over observations
        % same words --> difference '0'
        % diametral words --> sum 'wlen^wlen-1'
        T(jj) = sum((wseq(1,(startv+tt):(endv+tt)) - wseq(2,startv:endv)) == 0)/(endv - startv + 1);
        Tbar(jj) = sum((wseq(1,(startv+tt):(endv+tt)) + wseq(2,startv:endv)) == wlen^wlen - 1)/(endv - startv + 1);
    end
end

res = T - Tbar;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if pval > 0
    % significance
    if ensemble
        dtcrit = 2.7005*M^(-0.5179);
    else
        dtcrit = 2.7005*N^(-0.5179);
    end
    res(abs(res)<dtcrit) = 0;
    % quantiles
    res(abs(res)<quantile(res(:),1-pval)) = 0;
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% graphical output
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~silent
    if ensemble
        xres = max(abs(min(min((res)))),abs(max(max(res))));
        xres = xres + 0.1*xres + eps;
        figure;
        imagesc(res',[-xres xres]);
        colorbar;
        set(gca,'YTick',1:length(taurange));
        set(gca,'YTickLabel',taurange);
        title('pesct x -> y');
        grid on
    else
        figure;
        bar(res);
        set(gca,'XTick',1:length(taurange));
        set(gca,'XTickLabel',taurange);
        title('pesct x -> y');
        grid on
    end
end


