clear; close all; clc;

%% load data
load('data_1subj.mat');
[N,M]=size(So); 

%% parameters 
pfilter=0.94; %filter parameter for detrending
% MB parameters
q=20;
pmax=12; % maximum scanned model order
% MF parameters
m_max=10; % maximum embedding dimension
k=10; % number of neighbors
num_rnd=100; minshift=0; alpha_rnd=0.05; % CMI-based NUembed: surro-based termination criterion 

numsurr=100; % no of surrogates
perc=95; % confidence intervals calculation
minshift_surr = 20; % minimum shift for MF surrogate approach

pos_bar=[0.5 2 3.5]; pos_box=[1 2.5 4];

sel_targ{1}='AC'; sel_targ{2}='MAP'; sel_targ{3}='HP'; sel_targ{4}='RESP'; 
sel_source{1,1}='MAP'; sel_source{1,2}='HP'; sel_source{1,3}='RESP';
sel_source{2,1}='AC'; sel_source{2,2}='HP'; sel_source{2,3}='RESP';
sel_source{3,1}='AC'; sel_source{3,2}='MAP'; sel_source{3,3}='RESP';
sel_source{4,1}='AC'; sel_source{4,2}='MAP'; sel_source{4,3}='HP';

% ZeroLag, M*M matrix with 1 in the position ij when an instantaneous effect has been assigned from j to i
ZeroLag=[0 1 1 1; 0 0 1 1; 0 0 0 1; 0 0 0 0];

% init
EGC_dr_to_targ=nan(M,M-1);
EGC_dr_to_targ_s=nan(numsurr,M-1);
EGC_dr_to_targ_perc=nan(1,M-1);
sign=zeros(M,M-1);

%% pre-processing: AR filtering and removal of mean value
Sf=nan*ones(N,M); S=nan*ones(N,M);
for m=1:M
    Sf(:,m)=AR_filter(So,m,pfilter); % AR highpass filtered series
    S(:,m)=So(:,m)-mean(So(:,m)); % zero-mean series
end

%% MB analysis

% model order selection
pottaic = mos_idMVAR(S',pmax,0);
% MB extended conditional Granger causality
EGC=egc_gceMVAR_prior(S',pottaic,ZeroLag);
EGC=EGC./2; % corresponding information measure
for ii=1:M
    for jj=1:M
        if ii==jj; EGC(ii,jj)=nan; end
    end
end

figure('Windowstate','maximized')

for itarg=1:M
    EGC_dr_to_targ(itarg,:)=EGC(itarg,~isnan(EGC(itarg,:))); % vector 1x3 (no. of sources)

    % surrogate data analysis
    for ns=1:numsurr

        S_targ=surriaafft(S(:,itarg)); % random shuffling of the target
        Ss = S;
        Ss(:,itarg) = S_targ;

        % model order selection
        pottaic_s = mos_idMVAR(Ss',pmax,0);
        % MB extended conditional Granger causality
        EGCs=egc_gceMVAR_prior(Ss',pottaic_s,ZeroLag);
        EGCs=EGCs./2;
        for ii=1:M
            for jj=1:M
                if ii==jj; EGCs(ii,jj)=nan; end
            end
        end
        
        EGC_dr_to_targ_s(ns,:)=EGCs(itarg,~isnan(EGCs(itarg,:))); 

    end

    % test significance
    for isource=1:M-1

        % find the threshold
        EGC_dr_to_targ_perc(isource) = prctile(EGC_dr_to_targ_s(:,isource),perc);
        if EGC_dr_to_targ_perc(isource) < EGC_dr_to_targ(itarg,isource)
            sign(itarg,isource) = 1; 
        end

    end

    subplot(2,2,itarg)
    b=bar(pos_bar,EGC_dr_to_targ(itarg,:),0.3);
    ytips = b.YEndPoints;
    hold on;
    box=boxplot(EGC_dr_to_targ_s,'Positions',pos_box,'Widths',0.3);
    xticks(pos_bar); xticklabels(sel_source(itarg,:))
    ylabel('cGC [nats]');
    ylim([0 max(ytips)+0.1])
    for isource=1:M-1
        if sign(itarg,isource)==1
            text(pos_bar(isource),ytips(isource),'*',...
            'HorizontalAlignment','center','VerticalAlignment','bottom',...
            'Fontsize',18)
        end
    end
    title(horzcat('Target: ',sel_targ{itarg}))
end

sgtitle('Representative example: MB')

%% MF analysis

for m=1:M
    S_norm(:,m) = zscore(S(:,m));
end

cTE = nan(M);
sign_cTE = zeros(M);

figure('Windowstate','maximized')

for itarg = 1:M
    for isource = 1:M
    
        if isource ~= itarg
    
            V_tmp=m_max*ones(1,M);      
            candidates=its_SetLag(V_tmp,ones(1,M),ones(1,M),ZeroLag(itarg,:)); 
            out_nue=its_NUEknn(S_norm,itarg,candidates,k,num_rnd,minshift,alpha_rnd); 
            Vnu = out_nue.VL;
            out_cTE=its_PTEknn(S_norm,Vnu,isource,itarg,k);
            cTE(itarg,isource)=out_cTE.Txy_z;

            % surrogate data analysis
            for ns=1:numsurr
                maxshift=N-minshift_surr;
                lagshift=fix(rand(1,1).*(maxshift-minshift_surr+1)+minshift_surr);
                S_targ = circshift(S_norm(:,itarg),lagshift);        
                Ss = S_norm; Ss(:,itarg) = S_targ;
 
                out_nue_surr=its_NUEknn(Ss,itarg,candidates,k,num_rnd,minshift,alpha_rnd); 
                Vnu_surr = out_nue_surr.VL;
                out_cTE_surr=its_PTEknn(Ss,Vnu_surr,isource,itarg,k);
                cTE_surr(ns,isource)=out_cTE_surr.Txy_z;
            end

            % test significance
            th_cTE=prctile(cTE_surr(:,isource),perc);
            if cTE(itarg,isource) > th_cTE
                sign_cTE(itarg,isource) = 1;
            end

        end

    end

    cTE_tmp = cTE(itarg,:); cTE_tmp(itarg) = [];
    cTE_surr_tmp = cTE_surr; cTE_surr_tmp(:,itarg) = [];
    sign_cTE_tmp =  sign_cTE(itarg,:); sign_cTE_tmp(itarg) = [];

    subplot(2,2,itarg)
    b=bar(pos_bar,cTE_tmp,0.3);
    ytips = b.YEndPoints;
    hold on;
    box=boxplot(cTE_surr_tmp,'Positions',pos_box,'Widths',0.3);
    xticks(pos_bar); xticklabels(sel_source(itarg,:))
    ylabel('cTE [nats]');
    ylim([0 max(ytips)+0.1])
    for isource=1:M-1
        if sign_cTE_tmp(isource)==1
            text(pos_bar(isource),ytips(isource),'*',...
            'HorizontalAlignment','center','VerticalAlignment','bottom',...
            'Fontsize',18)
        end
    end
    title(horzcat('Target: ',sel_targ{itarg}))

end

sgtitle('Representative example: MF')