Description of vb_job_vb

0001 function vb_job_vb_dynamics(varargin)
0002 % Jointly estimate source activities and interactions
0003 %
0004 % [syntax]
0005 % vb_job_vb_dynamics(dbayes_parm)
0006 % vb_job_vb_dynamics(proj_root,dbayes_parm)    [old style]
0007 %
0008 % [input]
0009 % proj_root : <<string>> VBMEG project root directory.
0010 % dbayes_parm: <<struct>> Parameters for variational Bayes method.
0011 % --- fields of dbayes_parm
0012 %  (filename parameters) - relative path from proj_root (inside project files.)
0013 %  bayesfile  : <<string>> Current variance file used as initial values (.bayes.mat).
0014 %  dbayesfile : <<string>> File for saving source activities and interactions (.dbayes.mat).
0015 %
0016 %                        - fullpath is required (outside project files.)
0017 %  dmri_file  : <<string>> Diffusion MRI file (.dmri.mat).
0018 %
0019 %  (prior parameters: see Fukushima et al., 2015, Theory setion)
0020 %  ieta0 : <<double>> Prior mean of the inverse variance parameters of
0021 %          the MAR coefficients (\bar{\eta}_0^{-1} in TeX representation)
0022 %  g0    : <<double>> Shape parameter of gamma probability distribution
0023 %          for the inverse variance parameters of the MAR coefficients
0024 %
0025 %  (forward model parameters)
0026 %  area_key       : <<string>> ID of cortical area at which current
0027 %                   dipoles are assumed.
0028 %  area_key_global: <<string>> ID of cortical area at whith current
0029 %                   dipoles are assumed in estimation of baseline
0030 %
0031 %  (algorithm parameters)
0032 %  Ntrain    : <<int>> Iteration number of algorithm.
0033 %  Npre_train: <<int>> Iteration number of algorithm (=Ntrain).
0034 %
0035 % [output]
0036 % Estimation result will be saved to the file specified by
0037 % dbayes_parm.dbayesfile with the following variables:
0038 %
0039 %  Model.a    : current noise variance.
0040 %  Model.sx   : observation noise variance.
0041 %  Model.v    : set to 0.
0042 %  Model.Z    : current sources
0043 %  Model.Vz   : source covarinace matrices.
0044 %  Model.vz_sample : # of time samples for each cell entry of Model.Vz
0045 %  Model.mar  : MAR coefficients
0046 %  Model.MAR  : MAR matrix
0047 %  Model.MAR0 : MAR matrix at the 1st iterative update
0048 %  Model.Vmar : Covariace matrices of MAR coefficients
0049 %  Model.ieta : Variance parameters of MAR coefficients
0050 %  Model.a0   : prior mean of current noise variance.
0051 %  Model.sx0  : observation noise variance (baseline).
0052 %  Model.v0   : set to 0.
0053 %  Model.ix   : Dipole indices of Model.a.
0054 %  Model.Cov  : Noise covariance matrix.
0055 %  dbayes_parm: Settings specific to dynamic estimation.
0056 %  vb_parm    : Settings for VB iterative updates.
0057 %  Info       : Free energy and its constituents.
0058 %
0059 % [reference]
0060 % Fukushima et al. (2015) NeuroImage 105, 408-427.
0061 %
0062 % [history]
0063 %
0064 % 2015/03/18 M.Fukushima
0065 %
0066 % Copyright (C) 2011, ATR All Rights Reserved.
0067 % License : New BSD License(see VBMEG_LICENSE.txt)
0068 
0069 if length(varargin) == 1
0070   proj_root = [];
0071   dbayes_parm = varargin{1};
0072 elseif length(varargin) == 2
0073   proj_root = varargin{1};
0074   dbayes_parm = varargin{2};
0075 end
0076 
0077 %
0078 % --- Previous check
0079 %
0080 if ~isfield(dbayes_parm, 'bayes_file')
0081     error('dbayes_parm.bayes_file is a required field.');
0082 end
0083 init_value_file = fullfile(proj_root, dbayes_parm.bayes_file);
0084     
0085 if exist(init_value_file, 'file') ~= 2
0086     error('dbayes_parm.bayes_file not found.');
0087 end
0088 
0089 if ~isfield(dbayes_parm, 'dmri_file') || ...
0090    exist(dbayes_parm.dmri_file, 'file') ~= 2
0091     error('dbayes_parm.dmri_file not found.');
0092 end
0093 if ~isfield(dbayes_parm, 'dbayes_file') || isempty(dbayes_parm.dbayes_file)
0094     error('dbayes_parm.dbayes_file: Output filename not specified.');
0095 end
0096 if ~isfield('dbayes_parm', 'area_key')
0097     dbayes_parm.area_key = [];
0098 end
0099 if ~isfield('dbayes_parm', 'area_key_global');
0100     dbayes_parm.area_key_global = [];
0101 end
0102 
0103 if length(dbayes_parm.dbayes_file) < 11
0104     error('Please check output filename.(.dbayes.mat)');
0105 end
0106 save_file = fullfile(proj_root, dbayes_parm.dbayes_file);
0107 if ~strcmp(save_file(end-10:end), '.dbayes.mat')
0108     save_file = [save_file, '.dbayes.mat'];
0109     warning('save filename doesn''t have an extension(.dbayes.mat)\n filename was changed to %s', save_file);
0110 end
0111 
0112 %
0113 % --- Main Procedure
0114 %
0115 proj_root = vb_rm_trailing_slash(proj_root);
0116 
0117 const = vb_define_verbose;
0118 VERBOSE_LEVEL_INFO = const.VERBOSE_LEVEL_INFO;
0119 
0120 load(init_value_file, 'bayes_parm');
0121 if ~isfield(bayes_parm, 'spatial_smoother')
0122     dp = vb_set_bayes_default_parameters;
0123     bayes_parm.spatial_smoother = dp.spatial_smoother;
0124 end
0125 
0126 if ~strcmp(bayes_parm.forward_model,'focal')
0127     warning(['foward_model was changed from ' ...
0128              '''' bayes_parm.forware_model '''' ...
0129              ' to ''focal''.']);
0130     bayes_parm.forward_model = 'focal';
0131 end
0132 if bayes_parm.patch_norm == ON
0133     error(['The specified initial value file was calculated using ', ...
0134            'bayes_parm.patch_norm=ON.', sprintf('\n'), ...
0135            'Please calculate it using bayes_parm.patch_norm=OFF.']);
0136 end
0137 
0138 %
0139 % --- Overwrite settings
0140 %
0141 
0142 % bayes_parm.megfile
0143 % bayes_parm.megfile_baseline
0144 if isfield(dbayes_parm, 'meg_file') && ~isempty(dbayes_parm.meg_file)
0145    if ~ischar(dbayes_parm.meg_file)
0146        error('dbayes_parm.meg_file : accept only char string.');
0147    end
0148    bayes_parm.megfile          = {dbayes_parm.meg_file};
0149    bayes_parm.megfile_baseline = {dbayes_parm.meg_file};
0150 end
0151 
0152 %
0153 % Change path to absolute path
0154 %
0155 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0156 vb_disp('----- Dynamics estimation -----',VERBOSE_LEVEL_INFO)
0157 bayes_parm_abs = vb_parm_absolute_path(proj_root, bayes_parm);
0158 
0159 % create temp area file and add diffusion parcellation area
0160 % to the copied area file. then use it instead of original area file.
0161 area_key = 'dmri_parcellation_area';
0162 [bayes_parm_abs, dbayes_parm] = ...
0163     inner_replace_area_file_with_diffusion_area_file(...
0164                                     bayes_parm_abs, dbayes_parm, area_key);
0165 %
0166 % Observation noise estimation
0167 %
0168 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0169 vb_disp('----- Observation noise model specification', ...
0170         VERBOSE_LEVEL_INFO); 
0171 
0172 % Cov = Normalized noise covariance matrix
0173 [Cov,sx00] = vb_observation_noise_specification(bayes_parm_abs);
0174 
0175 %
0176 % Preparation of lead fields
0177 %
0178 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0179 vb_disp('----- Leadfield matrix preparation',VERBOSE_LEVEL_INFO); 
0180 
0181 lf_parm.brainfile  = bayes_parm_abs.brainfile;
0182 lf_parm.areafile   = bayes_parm_abs.areafile;
0183 lf_parm.patch_norm = bayes_parm_abs.patch_norm;
0184 lf_parm.expand_spatial_filter = bayes_parm_abs.expand_spatial_filter;
0185 lf_parm.spatial_smoother = bayes_parm_abs.spatial_smoother;
0186 
0187 % Global window
0188 vb_disp('--- Lead field matrix of global window ',VERBOSE_LEVEL_INFO);
0189 lf_parm.basisfile = bayes_parm_abs.basisfile_global;
0190 lf_parm.area_key  = dbayes_parm.area_key_global;
0191 lf_parm.reduce    = bayes_parm_abs.reduce_global;
0192 lf_parm.Rfilt     = bayes_parm_abs.Rfilt_global;
0193 
0194 [Gall, ix_all]    = vb_leadfield_preparation(lf_parm);
0195 
0196 % Focal window
0197 vb_disp('--- Lead field matrix of focal window ',VERBOSE_LEVEL_INFO);
0198 lf_parm.basisfile = bayes_parm_abs.basisfile;
0199 lf_parm.area_key  = dbayes_parm.area_key;
0200 lf_parm.reduce    = bayes_parm_abs.reduce;
0201 lf_parm.Rfilt     = bayes_parm_abs.Rfilt;
0202 lf_parm.remove_area_key = [];
0203 
0204 [Gact, ix_act] = vb_leadfield_preparation(lf_parm);
0205 
0206 % Note: Gall and Gact are not normalized here!
0207 
0208 %
0209 % Normalization of lead fields
0210 %
0211 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0212 Nsession = length(Gact);  % Number of sessions
0213 
0214 % Normalization constant of leadfield
0215 Bnorm = zeros(Nsession,1);
0216 for n = 1: Nsession
0217   GGall = sum(Gall{n}.^2, 2);
0218   Bnorm(n)  = sum(GGall)/length(GGall);
0219 end
0220 bsnorm = sqrt(mean(Bnorm));
0221 
0222 % Normalization of leadfield
0223 for n=1:Nsession
0224   nGall{n} = Gall{n}/bsnorm;
0225   nGact{n} = Gact{n}/bsnorm;
0226 end
0227 
0228 clear Gall Gact;
0229 
0230 %
0231 % MEG data preparation
0232 %
0233 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0234 [B, Ntrial, tmp, Tsample] = vb_megdata_preparation(bayes_parm_abs);
0235 
0236 %
0237 % Check variable consistency
0238 %
0239 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0240 vb_check_variable(B, Cov, nGact, nGall);
0241 
0242 %
0243 % Estimation of baseline current variance
0244 %
0245 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0246 vb_disp('----- Baseline current variance estimation',VERBOSE_LEVEL_INFO);
0247       
0248 % current baseline estimation at pretrigger period
0249 % together with spherical sensor noise estimation
0250 [v0, sx0] = vb_baseline_variance_acquisition(bayes_parm_abs, nGall);
0251 
0252 if strcmp(bayes_parm_abs.forward_model,'focal');
0253     sx0 = sx00;
0254 end
0255 
0256 %
0257 % Set VB parameters
0258 %
0259 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0260 
0261 Njact = size(nGact{1},2);
0262 Njall = size(nGall{1},2);
0263 vb_parm = vb_set_vb_parm(bayes_parm_abs, B, Njact, Njall, v0, sx0);
0264 vb_parm.Ntrain = dbayes_parm.Ntrain;
0265 
0266 %
0267 % Set prior parameters
0268 %
0269 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0270 prior_parm.ix = ix_act;
0271 prior_parm.v0 = v0;
0272 prior_parm.a0 = bayes_parm.a0;
0273 prior_parm.a0_act = bayes_parm.a0_act;
0274 
0275 w = bayes_parm.prior_weight;
0276 prior_parm.Ta0 = Ntrial * Tsample * w / (2*(1-w));
0277 prior_parm.Ta0_act = prior_parm.Ta0;
0278 
0279 prior_parm.act = vb_get_act(bayes_parm_abs.actfile, ...
0280                          bayes_parm_abs.act_key);
0281 [a0,Ta0] = vb_set_vb_prior(prior_parm);
0282 vb_parm.a0 = a0;
0283 vb_parm.Ta0 = Ta0;
0284 
0285 %
0286 % Wiener filter (removed from the previous version)
0287 %
0288 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0289 
0290 Model = [];
0291 
0292 %
0293 % Estimation with focal (and global) window
0294 %
0295 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0296 
0297 %
0298 % Structural Connectivity
0299 %
0300 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0301 
0302 connectivity_matrix = [];
0303 delay_ms_matrix     = [];
0304 load(dbayes_parm.dmri_file, 'connectivity_matrix', 'delay_ms_matrix');
0305 
0306 vb_parm.Ind = connectivity_matrix;
0307 meginfo = vb_load_meg_info(bayes_parm_abs.megfile{1});
0308 vb_parm.Delta = round(delay_ms_matrix*meginfo.SampleFreq*1e-3); % ms -> samples
0309 
0310 %
0311 % Hyperparameters for the variance parameter of the MAR coefficients
0312 %
0313 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0314 
0315 vb_parm.ieta0 = dbayes_parm.ieta0;
0316 vb_parm.g0 = dbayes_parm.g0;
0317 
0318 %  'focal',
0319 vb_parm.v0  = 0;
0320 vb_parm.Tv0 = 0; 
0321 
0322   vb_parm.a = vb_calc_initialvalue_a(init_value_file, ...
0323     bayes_parm_abs.brainfile, dbayes_parm.dmri_file, bsnorm, vb_parm.a_min);
0324 [Model,Info] = vbmeg_ard_estimate25_3(B, [], nGact, ...
0325                                    Cov, vb_parm, Model);
0326 
0327 Model.sx0 = sx0;
0328 
0329 %
0330 % Denormalization
0331 %
0332 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0333 
0334   %  Alpha is scaled back by bsnorm
0335   fprintf('Alpha is scaled back by bsnorm\n')
0336   Model.a = (Model.a)./bsnorm^2;
0337   Model.v = (Model.v)./bsnorm^2;
0338   Model.a0 = (vb_parm.a0)./bsnorm^2;
0339   Model.v0 = (vb_parm.v0)./bsnorm^2;
0340   Model.Cov = Cov;
0341   Model.ix = ix_act;
0342   for ns = 1:Nsession
0343     Model.Z{ns} = (Model.Z{ns})./bsnorm;
0344     for nd = 1:length(Model.Vz{ns})
0345       Model.Vz{ns}{nd} = (Model.Vz{ns}{nd})./bsnorm^2;
0346     end
0347   end
0348 
0349 %
0350 % Save result
0351 %
0352 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0353 
0354 fprintf('----- Save estimation result in %s\n',save_file)
0355 
0356 vb_fsave(save_file, ...
0357         'Model','dbayes_parm','vb_parm',...
0358          'Info');
0359 
0360 % % project_file save
0361 % proj_file = get_project_filename;
0362 % if isempty(proj_file)
0363 %     return;
0364 % end
0365 
0366 % project_file_mgr('load', proj_file);
0367 % project_file_mgr('add', 'bayes_parm', bayes_parm);
0368 
0369 function [bayes_parm_abs, dbayes_parm] = ...
0370     inner_replace_area_file_with_diffusion_area_file(...
0371                                                 bayes_parm_abs, ...
0372                                                 dbayes_parm, area_key)
0373 % copy specified area file as a temporary file and dMRI are into the file,
0374 % then return the filename.
0375 % [Input]
0376 %    bayes_parm_abs  : initial value is calculated by this bayes_parm.
0377 %    dnayes_parm     : dynamics bayes estimation parameter.
0378 %    area_key        : save dmri area into a temporary file by this name.
0379 % [Output]
0380 %    bayes_parm_abs.areafile : temporary area filename including dMRI area.
0381 
0382 if ~isfield(dbayes_parm, 'area_key') || isempty(dbayes_parm.area_key)
0383     dbayes_parm.area_key = area_key;
0384 end
0385 if ~isfield(dbayes_parm, 'area_key_global') || isempty(dbayes_parm.area_key_global)
0386     dbayes_parm.area_key_global = area_key;
0387 end
0388 
0389 % find area key
0390 found = false;
0391 load(bayes_parm_abs.areafile, 'Area');
0392 for k=1:length(Area)
0393     if strcmp(Area{k}.key, area_key)
0394         found = true;
0395     end
0396 end
0397 if found, return; end
0398 
0399 % bayes_parm_abs.areafile
0400 load(dbayes_parm.dmri_file, 'vbmeg_area_ix');
0401 
0402 % copy area file to temp area file and put dMRI area into the file.
0403 temp_area_file = [tempname(tempdir), '.area.mat'];
0404 copyfile(bayes_parm_abs.areafile, temp_area_file);
0405 area_new.Iextract = vbmeg_area_ix;
0406 area_new.key      = area_key;
0407 vb_add_area(temp_area_file, area_new);
0408 
0409 % use temp_area_file instead of original area file.
0410 bayes_parm_abs.areafile = temp_area_file;
vb_job_vb_dynamics

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
