forward_prediction

 Compute current sources predicted by a linear MAR model

 Input
     MAR: MAR matrix                 [Nvact x Nvact]
   Delta: Time lag matrix            [Nvact x Nvact]
   omega: Input signals to the model [Nvact x Nt]

 Output
   Jpred: Predicted current sources  [Nvact x Nt]

 2014/10/20 M.Fukushima

 Copyright (C) 2011, ATR All Rights Reserved.
 License : New BSD License(see VBMEG_LICENSE.txt)

0001 function [Jpred] = forward_prediction(MAR, Delta, omega)
0002 % Compute current sources predicted by a linear MAR model
0003 %
0004 % Input
0005 %     MAR: MAR matrix                 [Nvact x Nvact]
0006 %   Delta: Time lag matrix            [Nvact x Nvact]
0007 %   omega: Input signals to the model [Nvact x Nt]
0008 %
0009 % Output
0010 %   Jpred: Predicted current sources  [Nvact x Nt]
0011 %
0012 % 2014/10/20 M.Fukushima
0013 %
0014 % Copyright (C) 2011, ATR All Rights Reserved.
0015 % License : New BSD License(see VBMEG_LICENSE.txt)
0016 
0017 % # of sources
0018 Nvact = size(MAR,1);
0019 % # of time instances
0020 Nt = size(omega,2);
0021 % Binarized structural connectivity matrix
0022 Ind = (Delta~=0);
0023 
0024 % Connection numbers
0025 sum_ind = sum(Ind,2);
0026 % Time lags
0027 unidlt = full(setdiff(unique(Delta),0));
0028 
0029 % Store matrix entries into a smaller size of matrix
0030 Indx   = zeros(Nvact,max(sum_ind));
0031 Deltax = zeros(Nvact,max(sum_ind));
0032 MARx   = zeros(Nvact,max(sum_ind));
0033 for nv = 1:Nvact
0034   Indx(nv,1:sum_ind(nv))   = find(Ind(nv,:)==1);
0035   Deltax(nv,1:sum_ind(nv)) = Delta(nv,Delta(nv,:)~=0);
0036   MARx(nv,1:sum_ind(nv))   = MAR(nv,Ind(nv,:)~=0);
0037 end
0038 
0039 % Forward prediction by a linear MAR model
0040 Jpred = zeros(Nvact, Nt);
0041 for t = 1:Nt
0042   Jtmp = [zeros(Nvact, unidlt(end)) Jpred];
0043   t_end = t + unidlt(end);
0044   
0045   % Prediction [MEX ver.]
0046   AJpred = atimesj(MARx, Indx, Deltax, sum_ind, ...
0047     Jtmp(:,t_end:-1:(t_end-unidlt(end))));
0048 %   % Prediction [Matlab ver.]
0049 %   AJpred = zeros(Nvact,1);
0050 %   for nv = 1:Nvact
0051 %     for ii = 1:sum_ind(nv)
0052 %       AJpred(nv) = AJpred(nv) + ...
0053 %         MARx(nv,ii) * Jtmp(Indx(nv,ii),t_end-Deltax(nv,ii));
0054 %     end
0055 %   end
0056 
0057   % Predicted current sources
0058   Jpred(:,t) = AJpred + omega(:,t);
0059 end