Home > vbmeg > functions > estimation > bayes > vbmeg_current_reconstruct_tr.m

vbmeg_current_reconstruct_tr

PURPOSE ^

Current reconstruction for each trial using Bayesian inverse filter

SYNOPSIS ^

function [Jact ,ix_act, Jbck, ix_bck, Var, JactInfo] =vbmeg_current_reconstruct_tr(proj_root, bayesfile, ix_area,tsubsamp, jactdir, mode)

DESCRIPTION ^

 Current reconstruction for each trial using Bayesian inverse filter
 --- Syntax
 [Jact ,ix_act, Jbck, ix_bck, Varact, JactInfo] ...
  = vbmeg_current_reconstruct_tr(proj_root, bayesfile, ...
                                ix_area, tsubsamp, jactdir)
  This program output current files for each trial to specified directory
  If trial average current is needed, use 'vbmeg_current_reconstruct'

 --- Input
 bayesfile : result file obtained by VB estimation

 --- Optional Input
 ix_area : Position index to calculate estimated current
           If 'ix_area' is empty or not given, 
              currents in the active region are calculated
 tsubsamp : Specify subsampled time index
           If 'tsubsamp' is empty or not given, 
              time subsampling is not done
 jactdir  : output directory of current files for each trial
            relative path from proj_root
            do not include '/' or '\' in the end of string
            If 'jactdir' is not given,  jactdir = 'temp'

 ---Output
 Jact    : active current
 ix_act  : Vertex index corresponding to active current

 Jact(n,t,:) is the current at the vertex 'ix_act(n)' & the time 't'
 Jact(Nact,Nsample,Ntrial) 
   Nact     : # of active region, 
   Nsample  : # of time sample, 
   Ntrial   : # of trials 

 Jbck     : background current
 ix_bck   : Vertex index corresponding to background current
            same for 'Jact'
 Varact   : posterior variance of act_area
   When active area is empty,
   posterior variance for background area is calculated
 ----
 MEG data : B{Nsession}(Nsensors, Tsample, Ntrials)

 2004-9-25 M. Sato
 2004-12-5 M. Sato : Modified
 2004-12-22 M. Sato : Modified
 2005-1-5  M. Sato Modified
 2005-1-15  M. Sato Modified
 2005-08-19 O. Yamashita
 * Posterior variance is calculated simultaneously
 * Variables are renamed
 * Corresponding to new bayesfile format
 2005-08-19 O. Yamashita ver.30b
 2005-09-30 O. Yamashita
 * Minor bug fix when Lact >= 2 & Lact_var = 1
 2006-08-31 M. Sato
 * Remove focal area from global area
 * VBMEG result & liedfield is loaded by vb_load_result
 * Inverse filter calculation is done in vb_invfilter_calc
 2009-04-02 Taku Yoshioka
  Parameter name changed within this code for readability
  (just replacing 'resultfile' to bayesfile)

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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [Jact ,ix_act, Jbck, ix_bck, Var, JactInfo] = ...
0002     vbmeg_current_reconstruct_tr(proj_root, bayesfile, ix_area, ...
0003                   tsubsamp, jactdir, mode)
0004 % Current reconstruction for each trial using Bayesian inverse filter
0005 % --- Syntax
0006 % [Jact ,ix_act, Jbck, ix_bck, Varact, JactInfo] ...
0007 %  = vbmeg_current_reconstruct_tr(proj_root, bayesfile, ...
0008 %                                ix_area, tsubsamp, jactdir)
0009 %  This program output current files for each trial to specified directory
0010 %  If trial average current is needed, use 'vbmeg_current_reconstruct'
0011 %
0012 % --- Input
0013 % bayesfile : result file obtained by VB estimation
0014 %
0015 % --- Optional Input
0016 % ix_area : Position index to calculate estimated current
0017 %           If 'ix_area' is empty or not given,
0018 %              currents in the active region are calculated
0019 % tsubsamp : Specify subsampled time index
0020 %           If 'tsubsamp' is empty or not given,
0021 %              time subsampling is not done
0022 % jactdir  : output directory of current files for each trial
0023 %            relative path from proj_root
0024 %            do not include '/' or '\' in the end of string
0025 %            If 'jactdir' is not given,  jactdir = 'temp'
0026 %
0027 % ---Output
0028 % Jact    : active current
0029 % ix_act  : Vertex index corresponding to active current
0030 %
0031 % Jact(n,t,:) is the current at the vertex 'ix_act(n)' & the time 't'
0032 % Jact(Nact,Nsample,Ntrial)
0033 %   Nact     : # of active region,
0034 %   Nsample  : # of time sample,
0035 %   Ntrial   : # of trials
0036 %
0037 % Jbck     : background current
0038 % ix_bck   : Vertex index corresponding to background current
0039 %            same for 'Jact'
0040 % Varact   : posterior variance of act_area
0041 %   When active area is empty,
0042 %   posterior variance for background area is calculated
0043 % ----
0044 % MEG data : B{Nsession}(Nsensors, Tsample, Ntrials)
0045 %
0046 % 2004-9-25 M. Sato
0047 % 2004-12-5 M. Sato : Modified
0048 % 2004-12-22 M. Sato : Modified
0049 % 2005-1-5  M. Sato Modified
0050 % 2005-1-15  M. Sato Modified
0051 % 2005-08-19 O. Yamashita
0052 % * Posterior variance is calculated simultaneously
0053 % * Variables are renamed
0054 % * Corresponding to new bayesfile format
0055 % 2005-08-19 O. Yamashita ver.30b
0056 % 2005-09-30 O. Yamashita
0057 % * Minor bug fix when Lact >= 2 & Lact_var = 1
0058 % 2006-08-31 M. Sato
0059 % * Remove focal area from global area
0060 % * VBMEG result & liedfield is loaded by vb_load_result
0061 % * Inverse filter calculation is done in vb_invfilter_calc
0062 % 2009-04-02 Taku Yoshioka
0063 %  Parameter name changed within this code for readability
0064 %  (just replacing 'resultfile' to bayesfile)
0065 %
0066 % Copyright (C) 2011, ATR All Rights Reserved.
0067 % License : New BSD License(see VBMEG_LICENSE.txt)
0068 
0069 
0070 if ~exist('ix_area','var'),  ix_area = [];end;
0071 if ~exist('tsubsamp','var'), tsubsamp = [];end;
0072 if ~exist('jactdir','var'),  jactdir = 'temp'; end;
0073 if ~exist('mode','var'),     mode = 0; end;
0074 
0075 
0076 %%%%%% load result
0077 if ~isempty(proj_root)
0078     jactdir_ab = fullfile(proj_root, jactdir);
0079 else
0080     jactdir_ab = [jactdir];
0081 end  
0082 
0083 %check directory
0084 if exist(jactdir_ab) ~= 0
0085     msg = [ 'jact_dir : ' jactdir_ab ' exists. Do you over write?'];
0086     btn = questdlg(msg, 'confirm', 'Yes', 'No', 'Yes');
0087     if strcmp(btn, 'No')
0088         error('processing aborted.');
0089         return;
0090     end
0091 end
0092 
0093 if mkdir(proj_root, jactdir) ~= 1,  error('mkdir failed.'); end
0094 
0095 %%%%%% load result and MEG data
0096 [B, Model, nGact, nGall, nGGall, Cov, Wact, Wbck, ...
0097     bsnorm, ix_act, ix_bck, vb_parm] ...
0098     = vb_load_result(proj_root, bayesfile, ix_area);
0099 
0100 % B      : MEG data
0101 % Model  : Estimated model
0102 % nGact  : Normalized leadfield of focal window
0103 % nGall  : Normalized leadfield of global window
0104 % Cov    : Sensor noise covariance matrix
0105 % bsnorm : Normalizing constant
0106 % Wact   : Spatial smoothing matrix of focal window
0107 % Wbck   : Spatial smoothing matrix of global window
0108 % ix_act : Vertex index corresponding to active current
0109 % ix_bck : Vertex index corresponding to background current
0110 
0111 % # of vertex
0112 Nvact       = vb_parm.Nvact;         % # of active vertex
0113 Njact_area = vb_parm.Njact_area;
0114 Njall_area = vb_parm.Njall_area;
0115 
0116 %% Constant
0117 Nsession = vb_parm.Nsession;    % Number of sessions
0118 Ntrials  = zeros(Nsession,1);
0119 Tsample  = vb_parm.Tsample;        % Number of total time sample
0120 Twindow  = vb_parm.Twindow;
0121 
0122 % Temporal subsampling index
0123 if isempty(tsubsamp),
0124     tsubsamp = 1:Tsample;
0125 end
0126 
0127 % Temporal smoothing window weight
0128 [Tweight ,Tindex, Nindex, Nsample] = ...
0129     vb_overlapped_timewindow_weight(Twindow, Tsample, tsubsamp, mode);
0130 
0131 Nwindow = length(Nindex);       % # of time window
0132 
0133 % Nsample : # of subsampled data
0134 
0135 if mode == 1,fprintf('Non-overlapped concatenation mode\n'); end;
0136 
0137 %% Posterior variance
0138 if Njact_area > 0
0139     Var    = zeros(Njact_area,Nsample);
0140     Varact = zeros(Njact_area,Nsample);
0141 elseif Njall_area > 0
0142     Var    = zeros(Njall_area,Nsample);
0143     Varact = zeros(Njall_area,Nsample);
0144 end
0145 
0146 %% Current
0147 Jact_ave = zeros(Njact_area,Nsample);
0148 Jbck_ave = zeros(Njall_area,Nsample);
0149 
0150 %% Estimated current variance
0151 if isfield(Model,'a')
0152     a_inv = Model.a;
0153 else
0154     a_inv = zeros(Nvact,Nwindow);
0155 end
0156 if isfield(Model,'v')
0157     v_inv = Model.v;
0158 else
0159     v_inv = zeros(Nwindow,1);
0160 end
0161 
0162 %% Inverse filter
0163 KW  = cell(Nwindow,1);
0164 KW0 = cell(Nwindow,1);
0165 
0166 Ntrial_all = 0;
0167 
0168 % Session loop
0169 for n=1:Nsession
0170     Ntry = size(B{n},3);
0171     Ntrials(n) = Ntry;
0172     
0173     % Following notation is used in comments
0174     % G0 : original leadfield matrix (non-smoothed)   : []
0175     % A :  variance of active current                 : [A_inv]
0176     % tau : variance of active current                : [v_inv]
0177     % Wa : spatial smoothing matrix of active current : [Wact]
0178     % Wb : spatial smoothing matrix of background current : [Wback]
0179     % Ga : leadfield of active current (smoothed) G0*Wa   : [G]
0180     % Gb : leadfield of background current (smoothed) G0*Wb : []
0181     
0182     % Lead field for each session
0183     G     = nGact{n};   % Ga
0184     GbW  = nGall{n};   % Gb*Wb'
0185     GG   = nGGall{n};  % Gb*Gb'
0186     Covs = Cov{n};     % Sg
0187     Nch     = size(G,1);
0188     
0189     %%%% Calculate Inverse filter for current reconstruction
0190     for j=1:Nwindow
0191         % Time window loop
0192         % Subsampling time index
0193         Nid = Nindex{j};    % index in the total subsampled data
0194         Nt  = length(Nid);
0195         
0196         if isempty(Nid), continue; end;
0197 
0198         [KWt, KW0t, var_tmp] = vb_invfilter_calc(a_inv(:,j), v_inv(j), ...
0199                                   G, GbW, GG, Covs, Wact, Wbck);
0200         KW{j}  = KWt;
0201         KW0{j} = KW0t;
0202         var_tmp = Model.sx(j) * var_tmp./(bsnorm^2); 
0203         Varact(:,Nid) = Varact(:,Nid) + var_tmp * Tweight{j};
0204     end % Timewindow loop END
0205     
0206     %%%% Current reconstruction
0207     for m=1:Ntry
0208         Jact = zeros(Njact_area,Nsample);
0209         Jbck = zeros(Njall_area,Nsample);
0210         
0211         for j=1:Nwindow
0212             % Subsampling time index
0213             Tid = Tindex{j};    % subsampled time index
0214             Nid = Nindex{j};    % index in the total subsampled data
0215             
0216             if isempty(Nid), continue; end;
0217             
0218             weight = repmat(Tweight{j}, [Nch 1]);
0219             
0220             Bt = B{n}(:,Tid,m) .* weight; % Time window smoothing
0221             
0222             Jact(:,Nid) = Jact(:,Nid) + (KW{j} * Bt)/bsnorm;
0223             Jbck(:,Nid) = Jbck(:,Nid) + (KW0{j}* Bt)/bsnorm;
0224         end
0225         
0226         % Trial average current
0227         Jact_ave   = Jact_ave + Jact;
0228         Jbck_ave   = Jbck_ave + Jbck;
0229         Ntrial_all = Ntrial_all + 1;
0230         
0231         fname = sprintf('data_s%04dt%04d',n,m);
0232         vb_save([jactdir_ab '/' fname], 'Jact','Jbck', ...
0233                  'Jbck','ix_act','ix_bck','Varact');
0234         fprintf('progress... session:%04d/%04d, trial:%04d/%04d',...
0235                   n,Nsession,m,Ntry);
0236         %fprintf('%s',repmat(char(sprintf('\b')),1,46));
0237         fprintf('\n');
0238     end; % Trial loop END
0239     
0240     Var = Var + Ntry./max(Varact,eps);
0241 end; % Session loop END
0242 
0243 JactInfo.Ntrial   = Ntrials;
0244 JactInfo.Nsession = Nsession;
0245 JactInfo.jactdir  = jactdir;
0246 JactInfo.NJact    = Njact_area;
0247 JactInfo.Tsample  = Tsample;
0248 
0249 % Trial average current
0250 Jact = Jact_ave/Ntrial_all;
0251 Jbck = Jbck_ave/Ntrial_all;
0252 
0253 %Var  = Var/Ntrial_all;
0254 Var  = 1./Var;

Generated on Mon 22-May-2023 06:53:56 by m2html © 2005