%% EXAMPLE of MULTIVARIATE BLOCK ANALYSIS: realization of illustrative example and comparison between our multivariate measures and Geweke measures
% L.Faes and G. Nollo, "Measuring Frequency Domain Granger Causality for Multiple Blocks of Interacting Time Series", Biological Cybernetics 2013. DOI: 10.1007/s00422-013-0547-5
% (this code realizes Fig.7 of the paper)
% L.Faes, S. Erla and G. Nollo, "Block Partial Directed Coherence: a New Tool for the Structural Analysis of Brain Networks", International Journal of Bioelectromagnetism, Vol. 14, No. 4, pp. 162 - 166, 2012


clear; close all; clc;
nfft=512; % number of frequency bins
fc=1; % sample frequency
Q=8; % number of series
Mv=[2 2 3 1]'; % number of series in each block
pmax=2; % maximum lag

kiv=[1 2; 1 3; 3 4]';  % tested directions: 1->2, 1->3, 3->4
numsimu=4; % paper: 100 realizations
N=500; %length of simulated time series


%% MVAR coefficients
ca=[0.5 0.5];
cb=[0.4 0.4];
cc=[0.4 0.4 0.4];
cd=[0.7];
ck=[-0.5 0.5 0.5 0.5];

pmax=2; % maximum lag
r1=0.9; f1=0.1; % oscillation channel 1
r3=0.9; f3=0.25; % oscillation channel 3

a1=ca(1); a2=ca(2);
b1=cb(1); b2=cb(2);
c1=cc(1); c2=cc(2); c3=cc(3);
d1=cd(1);
k1=ck(1); k2=ck(2); k3=ck(3); k4=ck(4);
Su=eye(Q); %Residual covariance matrix (DIAGONAL)
B0=zeros(Q,Q); %Matrix of instantaneous effects:no instantaneous effects
Bk=NaN*ones(Q,Q,pmax);  
%effects at lag 1
Bk(1,:,1)=[2*r1*cos(2*pi*f1) 0 0 0 0 0 0 0];
Bk(2,:,1)=[k1 0 0 0 0 0 0 0];
Bk(3,:,1)=[0 0 2*r3*cos(2*pi*f3) 0 0 0 0 d1];
Bk(4,:,1)=[0 a2 0 0 0 0 0 0];
Bk(5,:,1)=[0 0 c1 0 0 0 0 0];
Bk(6,:,1)=[0 0 0 0 k2 0 0 0];
Bk(7,:,1)=[0 0 c3 0 0 k4 0 0];
Bk(8,:,1)=[0 0 0 b2 0 0 0 0];
%effects at lag 2
Bk(1,:,2)=[-r1^2 0 0 0 0 0 0 0];
Bk(2,:,2)=[0 0 0 0 0 0 0 0];
Bk(3,:,2)=[0 a1 -r3^2 0 0 0 0 0];
Bk(4,:,2)=[0 0 0 0 0 0 0 0];
Bk(5,:,2)=[0 0 0 0 0 k3 0 0];
Bk(6,:,2)=[0 0 c2 0 0 0 0 0];
Bk(7,:,2)=[0 0 0 0 0 0 0 0];
Bk(8,:,2)=[0 0 b1 0 0 0 0 0];
% concateno in matrice Bm
Bm=[];
for kk=1:pmax
    Bm=[Bm Bk(:,:,kk)];
end
Am=Bm; % no inst eff
E=eye(Q*pmax);AA=[Am;E(1:end-Q,:)];lambda=eig(AA);
lambdamax=max(abs(lambda));


%% Theoretical spectral functions
M=length(Mv);
Q=size(Am,1);
p=size(Am,2)/Q;
%%% theoretical block functions
[bDC,bPDC,mF,mG,bS,bP,bH,bAf,f] = block_fdMVAR(Am,Su,Mv,nfft,fc);

%% realization of the process
emFn=NaN*ones(nfft,numsimu,size(kiv,2)); emGn=emFn;
gewFn=emFn; gewFcn=emFn;
for cnt=1:numsimu
    disp(['simulation ' int2str(cnt) ' of ' int2str(numsimu)]);

    U=randn(Q,N); % uncorrelated gaussian innovations
    [Y]=MVARfilter(Am,U); % function from eMVAR toolbox

    % model identification
    [eAm,eSu,Yp,Up]=idMVAR(Y,p,0); % function from eMVAR toolbox
    %%% estimated block functions
    [ebDC,ebPDC,emF,emG,ebS,ebP,ebH,ebAf,f] = block_fdMVAR(eAm,eSu,Mv,nfft,fc);
   
    for k=1:size(kiv,2)
        emFn(:,cnt,k)=squeeze(emF(kiv(2,k),kiv(1,k),:));
        emGn(:,cnt,k)=squeeze(emG(kiv(2,k),kiv(1,k),:));
        [gewF,gewFc,f] = Geweke_f(Y,Mv,p,kiv(1,k),kiv(2,k),nfft,fc);
        gewFn(:,cnt,k)=gewF;
        gewFcn(:,cnt,k)=gewFc;
    end
    
end
emFn_m=mean(emFn,2); emFn_sd=std(emFn,0,2);
emGn_m=mean(emGn,2); emGn_sd=std(emGn,0,2);
gewFn_m=mean(gewFn,2); gewFn_sd=std(gewFn,0,2);
gewFcn_m=mean(gewFcn,2); gewFcn_sd=std(gewFcn,0,2);

%% graphs & disp
figure(1);
for k=1:size(kiv,2)
    subplot(size(kiv,2),2,k*2-1);
    plot(f,emFn_m(:,1,k),'r'); hold on; plot(f,emFn_m(:,1,k)+emFn_sd(:,1,k),'r:'); plot(f,emFn_m(:,1,k)-emFn_sd(:,1,k),'r:');
    axis([0 fc/2 -0.1 3]); title(['F^(^m^)^  ' int2str(kiv(1,k)) '->' int2str(kiv(2,k))]);
    
    subplot(size(kiv,2),2,k*2);
    plot(f,emGn_m(:,1,k),'r'); hold on; plot(f,emGn_m(:,1,k)+emGn_sd(:,1,k),'r:'); plot(f,emGn_m(:,1,k)-emGn_sd(:,1,k),'r:');
    axis([0 fc/2 -0.1 3]); title(['G^(^m^)^ ' int2str(kiv(1,k)) '->' int2str(kiv(2,k))]);    
end

figure(2);
for k=1:size(kiv,2)
    subplot(size(kiv,2),2,k*2-1);
    plot(f,gewFn_m(:,1,k),'k'); hold on; plot(f,gewFn_m(:,1,k)+gewFn_sd(:,1,k),'k:'); plot(f,gewFn_m(:,1,k)-gewFn_sd(:,1,k),'k:');
    axis([0 fc/2 -0.1 3]); title(['_g_e_w_ F  ' int2str(kiv(1,k)) '->' int2str(kiv(2,k))]);
    
    subplot(size(kiv,2),2,k*2);
    plot(f,gewFcn_m(:,1,k),'k'); hold on; plot(f,gewFcn_m(:,1,k)+gewFcn_sd(:,1,k),'k:'); plot(f,gewFcn_m(:,1,k)-gewFcn_sd(:,1,k),'k:');
    axis([0 fc/2 -0.1 3]); title(['_g_e_w_ Fc ' int2str(kiv(1,k)) '->' int2str(kiv(2,k))]);    
end
